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SECTION  1 


INTRODUCTION 

Treating  the  transient  dynamic  interaction  between  a  structure  in  contact  with  a  fluid  or 
elastic  medium  is  a  formidable  task.  Given  the  dynamical  equations  for  the  structure  and  a 
specification  of  the  initial  conditions,  external  dynamic  forces  and/or  incident-wave  field,  a 
doubly  asymptotic  approximation  (DAA)  provides  a  link  that  greatly  simplifies  the  analysis.  This 
simplification  allows  the  analyst  to  devote  most  of  his/her  computational  resources  to  the 
structural  model,  which  is  the  focus  of  interest,  by  minimizing  the  resources  required  for 
modelling  the  medium,  which  is  rarely  of  interest 

1.1  MOTIVATION 

In  the  1970's,  DAA's  were  first  developed  to  treat  the  acoustic  fluid-structure  interaction 
in  underwater  shock  problems  (Geers,  1971,  1974,  1978).  These  approximations  approach 
exactness  in  both  the  early-time/high-frequency  and  late-time/low-frequency  limits;  hence  the 
name  doubly  asymptotic.  Acoustic  DAA's  have  been  incorporated  in  a  variety  of  production 
computer  programs  that  are  routinely  used  for  engineering  analysis  (Ranlet,  et  al.,  1977; 
DeRuntz,  etal.,  1980;  DeRuntz  and  Brogan,  1980;  Neilson,  etal.,  1981;  Vasudevan  and  Ranlet, 
1982;  Atkatsh,  et  al.,  1987).  In  the  1980's,  the  acoustic  DAA  methodology  was  improved 
(Felippa,  1980;  Geers  and  Felippa,  1983;  Nicolas-Vullierme,  1989)  and  the  DAA  concept  was 
extended  to  elastodynamics  and  electromagnetics  (Underwood  and  Geers,  1981;  Mathews  and 
Geers,  1987;  Geers  and  Zhang,  1988). 

In  solids  and  fluids,  the  general  approach  has  been  to  regard  the  stress  and  displacement 
fields  in  the  medium  as  the  sum  of  those  associated  with  the  incoming  incident  wave  (if  there  is 
one)  and  those  associated  with  the  outgoing  scattered  wave  (of  which  there  is  always  one--if  there 
is  no  incident  wave,  the  scattered  wave  is  usually  called  the  radiated  wave).  Compatibility  of 
surface  tractions  and  displacements  provides  all  of  the  remaining  relations  needed  save  one: 
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a  relation  between  the  scattered-wave  stress  and  the  scattered-wave  displacement  over  the  surface 
of  the  structure  in  contact  with  the  medium. 

The  Kirchhoff  retarded-potential  integral  for  an  acoustic  medium  and  the  dynamic 
Somigliana  identity  for  an  elastic  medium  provide  exact  relations  connecting  scattered-wave 
tractions  and  displacements.  Unfortunately,  these  relations  are  integral  equations  over  the  contact 
surface  that  involve  field  variables  with  retarded-time  arguments;  hence  they  are  local  neither  in 
space  nor  in  time,  and  they  are  complicated.  These  characteristics  mitigate  against  computational 
efficiency,  prompting  the  development  of  simpler  relations.  Singly  asymptotic  approximations 
have  been  developed  that  apply  either  at  early  time  or  at  late  time  (but  not  both),  but  they  are 
not  sufficiently  robust  for  diverse  application.  In  contrast,  doubly  asymptotic  approximations  for 
external  domains  have  been  found  easy  to  use  and  remarkably  accurate  in  a  broad  apectrum  of 
applications. 

Recently,  interest  has  developed  in  DAA’s  for  internal  acoustic  domains,  motivated  by 
the  following  factors:  (1)  Internal  domains  of  practical  interest  often  possess  exceedingly  complex 
geometries,  which  makes  3-D  mesh  generation  for  finite-element  modelling  costly  and 
cumbersome;  (2)  Numerical  simulations  of  discontinuous  wave  fronts  through  3-D  finite-element 
meshes  are  typically  plagued  by  non-physical  osicillations,  which  compromise  the  value  of  the 
calculations:  (3)  An  exact  boundary-element  treament  based  on  Kirchhoffs  retarded  potential 
formulation  would  be  computationally  intensive,  typically  usurping  resources  that  are  needed  for 
accurate  structural  modelling. 

At  the  same  time,  interest  continues  in  the  advancement  of  DA  A  technology  for  ground 
shock  analysis;  of  particular  interest  is  the  extension  of  the  existing  first-order  DAA  for  the 
infinite  elastic  medium  to  treat  the  semi-infinite  elastic  medium.  Because  the  constitutive 
behavior  of  soil  is  so  often  highly  nonlinear,  the  placement  of  a  DAA  boundary  directly  on  the 
soil-structure  contact  surface  is  not  advisable.  However,  the  use  of  such  a  boundary  as  a  non¬ 
reflecting  boundary  at  a  modest  distance  from  the  contact  surface  is  most  attractive.  As  discussed 
by  Mathews  and  Geers,  1987,  a  DAA  nonreflecting  boundary  is  superior  to  the  singly  asymptotic 
boundaries  currently  used  in  most  codes. 

This  report  documents  recent  advances  in  DAA  technology.  First,  the  methodology  of 
formulating  DAA’s  is  systematized,  which  is  essential  for  the  development  of  high-order 
approximations.  Second,  first-  and  second-order  DAA’s  are  formulated  for  an  internal  acoustic 
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medium.  Third,  the  internal  DAA’s  are  evaluated  by  comparing  DAA  solutions  with  exact 
solutions  for  a  canonical  underwater-shock  problem;  because  solutions  to  this  problem  did  not 
previously  exist,  the  exact  solutions  are  provided  herein.  Fourth,  a  first-order  DAA  for  an  elastic 
half-space  is  formulated.  Fifth,  this  DAA  is  implemented  in  a  boundary-element  code  and 
numerical  results  for  two  canonical  problems  are  compared  with  results  currently  in  the  literature. 

1.2  REPORT  OUTLINE. 

Section  2  of  this  report  contains  a  review  of  the  first-order  DAA  (DAA,)  for  an  external 
acoustic  medium  featuring  the  method  of  operator  matching.  Both  integral-operator  and  matrix 
formulations  arc  presented,  and  a  modal  analysis  of  the  two  formulations  is  performed. 

DAA,  for  an  internal  acoustic  domain  is  formulated  in  Section  3.  The  separation  of  low- 
frequency  fluid  motion  into  dilatational  and  equivoluminal  components  is  shown  to  be  essential 
to  the  formulation.  Operator,  matrix  and  modal  developments  are  given,  all  based  on  operator 
matching. 

Section  4  contains  a  straightforward  formulation  of  the  second-order  DAA  (DAAj)  for  an 
external  acoustic  medium  produced  by  operator  matching.  This  extends  the  work  of  Felippa 
(1980a)  and  Nicolas-Vullieime  (1989),  avoiding  the  introduction  of  an  impedance  formalism  and 
retaining  the  advantages  of  Laplace  transformation.  The  corresponding  matrix  formulation  is  also 
presented,  but  a  modal  analysis  is  not,  because  the  matched  DAA2  does  not  diagonalize. 

The  matched  DAA2  for  an  internal  acoustic  domain  is  formulated  in  Section  5.  Again, 
the  separation  of  low-frequency  motion  into  dilatational  and  equivoluminal  components  is  central. 
Both  operator  and  matrix  forms  are  given,  but  uncoupled  modal  analysis  is  not  admissable. 

Section  6  describes  the  specialization  of  the  four  matched  DAA’s  to  axisymmetric  flow 
outside  and  inside  a  spherical  surface .  For  this  classical  geometry,  even  the  second-order  DAA’s 
submit  to  uncoupled  modal  decomposition  in  terms  of  Legendre  polynomials.  This  yields  modal 
DAA  equations  for  each  generalized  harmonic.  Also  provided  in  this  section  arc  exact  modal 
equations,  which  are  substantially  more  complicated  than  the  DAA  equations. 

In  Section  7.  exact  modal  response  equations  are  derived  for  a  previously  unsolved 
canonical  problem,  namely,  the  response  of  a  fluid-filled,  submerged  spherical  shell  to  a  transient 
acoustic  wave.  The  modal  equations  are  formulated  by  the  residual  potential  method  (Geers, 
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1969,  1971,  1972)  and  are  solved  by  numerical  integration  in  time.  Physical  responses  are  then 
obtained  by  modal  superposition.  Difficulties  with  poor  modal  convergence  are  successfully 
treated  by  obtaining  partial  closed-form  solutions  and  using  the  Cesaro  sum  (Apostol,  1957).  The 
numerical  solutions  thus  obtained  serve  as  basis  for  evaluating  the  internal  DAA’s  developed  in 
Sections  3  and  5. 

Numerical  results  for  the  fluid-filled,  submerged  spherical  shell  excited  by  a  plane  step- 
wave  are  presented  in  Section  8.  Exact,  DAAt,  and  DAA2  results  are  compared  to  assess  the 
accuracy  of  the  internal  DAA’s.  Also,  the  shock  response  of  the  fluid-filled  shell  is  contrasted 
with  that  of  an  empty  shell. 

In  Section  9.  systematic  DAA,  formulations  are  given  for  infinite  and  semi-infinite  elastic 
media,  both  in  operator  and  matrix  form.  Implementation  in  a  boundary-element  code  is 
described,  and  numerical  results  for  two  canonical  problems  are  compared  with  corresponding 
results  in  the  literature. 

Section  10  concludes  the  report  by  summarizing  the  work  conducted  and  listing  the 
principal  conclusions  reached  during  the  study. 

1.3  TECHNOLOGY  TRANSFER. 

The  implementation  of  external  acoustic  DAA’s  in  production  shock-analysis  codes  has 
improved  the  engineering  design  and  analysis  of  many  naval  structures.  The  internal  acoustic 
DAA2  formulated  in  Section  5  is  shown  in  Section  8  to  be  sufficiently  accurate  to  warrant  its 
early  implementation  in  those  codes.  In  the  meantime,  it  is  appropriate  to  seek  improved  internal 
DAA’s  in  order  to  raise  the  level  of  accuracy  to  that  exhibited  by  the  external  DAA2. 

More  research  is  needed  before  elastic  DAA’s  are  ready  for  production  analysis.  The 
first-order  DAA’s  for  infinite  and  semi-infinite  half-spaces  are  only  marginally  accurate,  which 
calls  for  the  development  of  second-order  DAA’s.  Fortunately,  such  development  can  make  good 
use  of  the  formulation  techniques  used  to  develop  acoustic  DAA’s. 

DAA’s  can  be  formulated  for  shock  response  analyses  involving  other  media,  such  as 
layered  media,  porous  media,  and  air  at  moderate  pressures;  higher-order  DAA’s  for 
electromagnetic  scattering  also  hold  promise.  What  was  orignally  developed  as  a  method  focused 
on  underwater  shock  analysis  is  emerging  as  one  of  substantially  broader  scope. 
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SECTION  2 


DAAt  FOR  AN  EXTERNAL  ACOUSTIC  DOMAIN 

Although  the  first-order  doubly  asymptotic  approximation  for  an  external  acoustic 
domain  was  given  some  twenty  years  ago  (Geers.  1971),  a  review  is  appropriate  here,  for  two 
reasons.  First,  such  a  review  provides  the  clearest  picture  of  the  DAA  concept,  and  second, 
it  introduces  the  operator  matching  method  at  the  simplest  level. 

2.1  RETARDED  POTENTIAL  FORMULATION. 

With  the  acoustic  pressure  p(r,t)  and  fluid-particle  displacement  u(r.t)  given  in  terms 
of  a  velocity  potential  0(r,t)  as 


p(r,t)  -  p0(r.t) 

(2.1) 

u(r.t)  -  -V0(r,t) 

where  p  is  the  mass  density  of  the  fluid,  an  overdot  denotes  differentiation  in  time,  and  V  is 
the  gradient  operator,  the  wave  equation  for  a  uniform  acoustic  fluid  is  (see,  e.g.. 
Pierce,  1981) 


c2V20  -  0  (2.2) 

where  c  is  the  speed  of  sound  in  the  fluid  and  V2  is  the  Laplacian  operator. 

With  n  as  the  normal  going  into  the  fluid  at  a  point  on  a  surface  S  that  bounds  the  fluid 
domain,  the  inward  fluid-particle  displacement  normal  to  that  surface  is  defined  by 
u  -  u-n. 

An  exact,  integral-equation  solution  to  (2.2)  is  given  by  Kirchhoff's  retarded  potential 
formulation  (RPF)  (see,  e.g..  Baker  and  Copson,  1939,  and  Sobolev,  1964),  which  may  be 
written  for  points  P  and  Q  on  S 


5 


2rr pp(t)  -  |  {pRpQUg(tR)  -  RPq  cos^RhIPq^r)  +  c'^pqPq^r)]}  *3Sq 

Js 


(2.3) 


where  RPq  -  |  rP  -  rg| .  <f>Rn  is  the  angle  between  RPq  and  n.  and  tR  is  the  retarded  time  t  - 
c-'Rpq;  the  line  through  the  integral  sign  indicates  that  the  point  P  is  excluded  from  the 
integral.  The  constant  2n  multiplying  pP(t)  on  the  left  side  of  (2.3)  indicates  that  the  point  P 
is  located  on  a  smooth  portion  of  the  surface  S.  If  P  is  not  on  the  surface,  but  is  inside  (or 
outside)  the  fluid  domain,  the  multiplying  constant  becomes  4 w  (or  zero);  if  S  is  not  smooth  at 
P.  but  instead  has  an  edge  or  a  corner  there,  the  multiplying  constant  becomes  the  value  of 
the  solid  angle  subtended  by  the  edge  or  corner. 

For  our  purposes,  it  is  convenient  formally  to  incorporate  the  singular  contribution 
2?rpP(t)  into  the  spatial  integral  of  (2.3)  and  then  take  the  Laplace  transform  of  the  result  to 
obtain 


Rpq  cos <ARn ( 1  +Rpq s/ c)  e^PQ5^  pQ (s)  dSQ  -  p  RP'Q  s2  e^PQ^u^s)  dSQ  (2.4) 

Js  Js 

2.2  FIRST-ORDER  EARLY-TIME  APPROXIMATION;  ETA,. 

Early-time  approximations  are  the  inverse  Laplace  transforms  of  algebraic  equations  to 
which  (2.4)  reduces  when  Rmaxs/c  »  1,  which  corresponds  in  the  time  domain  to  t  « 
c_,Rmax  (Geers,  1975).  The  first-order  ETA  for  an  external  acoustic  medium  was  first 
utilized  by  Mindlin  and  Bleich,  1953,  for  a  problem  in  polar  coordinates.  In  a  systematic 
analysis,  Felippa,  1980,  found  the  same  result  for  a  general  smooth  surface,  which  is,  in 
transform  space, 

ETA,(s):  pP(s)  -  pcsuP(s)  (2.5) 

Inverse  Laplace  transformation  yields  as  the  first-order  ETA  in  the  time  domain 

ETA,(t);  pP(t)  -  pcuP(t)  (2.6) 
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ETA,  is  clearly  a  local  approximation  in  space,  stating  that  each  element  of  the  surface 
S  independently  generates  a  plane  wave  that  propagates  normally  into  the  fluid.  In  addition, 
it  applies  equally  well  to  either  an  external  or  an  internal  fluid  domain.  Finally,  because  it 
approaches  exactness  only  as  s  -*  oo,  it  is  singly  asymptotic.  In  the  literature,  ETA,  is  often 
referred  to  as  the  plane  wave  approximation. 

2.3  FIRST-ORDER  LATE-TIME  APPROXIMATION:  LTA,. 

Late-time  approximations  are  the  inverse  Laplace  transforms  of  integral  equations  in 
space  to  which  (2.4)  reduces  when  R^^s/c  «  1.  which  corresponds  in  the  time  domain  to  t 
»  c-,Rmax  (Geers.  1975).  The  first-order  LTA  for  an  external  acoustic  domain  was  first 
utilized  by  Chertock,  1972.  It  may  be  readily  obtained  by  merely  expanding  the  exponentials 
in  (2.4)  in  a  Taylor's  series  as 

e-RpQs/c  _  ,  .  RpgS/c  +  i  (RpQs/c)2  -  ...  (2.7) 

Introducing  this  into  (2.4)  and  keeping  only  terms  of  order  s°  on  the  left  and  s2  on  the  right, 
we  obtain  LTA,  in  transform  space 


cosifon  pQ(s)  dSQ 


1  R' 


-l  -2 

PQ  S 


uQ(s)  dSQ 


(2.8) 


Note  that,  unlike  ETA,,  LTA,  is  not  spatially  local,  and  that,  like  ETA,,  it  is  singly 
asymptotic,  but  in  the  limit  s  -*  0  instead  of  s  ■*  oo.  In  the  literature,  LTA,  is  often  referred 
to  as  the  added  mass  approximation  or  the  virtual  mass  approximation. 

With  the  spatial  operator  definitions 


(2.9) 
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7Qq  s  Rpq  cos  0r„  qQ  dSQ 

Js 


LTA,  for  an  external  domain  may  be  expressed  in  transform  space  as 

LTA,(s):  /S~17Pq(s)  -  ps*up(s)  (2.10) 

or  in  the  time  domain  as 

LTA,(t):  r'7PQ(t)  -  pup(t)  (2.11) 

Here,  p denotes  the  inverse  of  the  operator  /?,  i.e.,  if  qP  produces  pP  «  /3qg  through  the  first 
of  (2.9),  then  pP  produces  qP  through  the  relation  qP  -  P~1Pq.  It  can  be  shown  that  0  is 
invertible. 

By  taking  the  Laplace  transform  of  (2.2),  considering  V2  ~  R^ax,  and  then  letting 
Rmaxs/c  -*  0,  one  readily  deduces  that  LTA,  constitutes  the  integral-equation  solution  to 
Laplaces  equation,  V20  -  0.  This  means  that  LTA,  pertains  to  the  irrotational  flow  of  an 
inviscid,  incompressible  fluid. 

2.4  FIRST-ORDER  DOUBLY  ASYMPTOTIC  APPROXIMATION:  DAA,. 

An  approximation  that  naturally  reduces  to  ETA,,  (2.5),  at  early  times  (s  -*  oo)  and  to 
LTA,,  (2.10),  at  late  times  (s  -  0)  is 

CsJ+2pp(s)  +  7Pq(s)  -  CpcsJ+3uP(s)  +  psJ/Suq(s)  (2.12) 

where  C  is  an  arbitrary  constant  and  j  ^  0.  This  approximation  has  two  flaws:  the  constants 
C  and  j  are  undetermined  and  the  inverse  transform  would  possess  derivatives  higher  than 
necessary.  Hence  we  reject  it  as  a  first-order  DAA. 

An  examination  of  (2.5)  and  (2.10)  reveals  that  a  relation  with  one  term  in  s°p(s)  and 
another  in  s‘p(s)  on  the  left,  and  with  one  term  in  s*u(s)  on  the  right,  is  capable  of  reducing  to 
the  two  singly  asymptotic  relations  in  the  appropriate  limits.  Hence,  as  the  first  step  in  the 
method  of  operator  matching ,  we  introduce  the  DAA,  trial  equation 
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[sP,  +  CP0]  Pq(s)  -  pC S* Up (s) 


(2.13) 


where  P0  and  P,  are  spatial  operators  (not  functions  of  s!).  For  s  -*  0,  we  write  this  equation 
as 


[/»„  +  0(s)]pq(s)  -  psJuP(s)  (2.14) 

and  match  it  to  (2.10)  as  s  -*■  0.  which  yields  P0  -  For  s  ■*  oo,  we  divide  (2.13)  through 

by  s,  write  the  result  as 


[/*,  +  0<S-1)]Pq(s)  -  pC  S  Up  (s)  (2.15) 

and  match  it  to  (2.5)  as  s  ■*  oo,  which  yields  F1pQ(s)  -  pp(s).  The  introduction  of  these  results 
into  (2.13)  produces  the  first-order  DAA  for  an  external  acoustic  domain,  expressed  in 
transform  space  as 

DAA,(s):  s pp (s)  +  c/r1ypQ(s)  -  pcs*uP(s)  (2.16) 

and  in  the  time  domain  as 


DAA,(t):  pp(t)  +  c/S~1ypQ(t)  -  pcuP(t)  (2.17) 

We  note  that,  as  might  be  anticipated.  DAA!  is  not  a  spatially  local  approximation. 


2.5  MATRIX  DAA,  FOR  BOUNDARY  ELEMENT  ANALYSIS. 


The  boundary  element  method  has  become  a  powerful  tool  for  obtaining  solutions  to 
problems  involving  complex  geometries  (see.  e.g..  Banner jee,  1981).  The  method  may  be 
described  with  considerable  generality  as  Petrov-Galerkin  finite-element  discretization  over 
the  boundary  of  a  spatial  domain  (Hughes.  1986).  To  use  the  method,  we  first  discretize  the 
pressure  and  normal-displacement  fields  on  the  surface  S  as 


Pq  (1)  -  vj  p(t) 


(2.18) 
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UQ(t)  -  vj  u(t) 


where  Vq  is  the  column  vector  of  shape-functions,  the  superscript  T  denotes  vector 
transposition,  and  p(t)  and  u(t)  are,  respectively,  the  column  vectors  for  nodal-pressure  and 
nodal-displacement  response.  To  be  able  to  represent  a  constant  field,  we  require  that  VqI  - 
1,  where  1  is  the  unit  vector. 

Next,  we  "preoperate"  the  DAAi  equation,  (2.17),  through  by  the  operator  /3.  insert 
(2.18),  and,  with  a  column  vector  of  weight-functions  wP,  form  the  weighted-residual 
equations 


|  wp  dSP  p(t)  +  c  |  wP  yv£  dSP  p(t)  -  pc  J  wP  0v^  dSP  u(t)  (2.19) 


which  can  be  written  more  compactly  as 

Bp(t)  +  cCp(t)  -  pcBu(t)  (2.20) 


where,  from  (2.9), 


C 


■  IslsWp  R 


P  ^PQ  VQ  dSp 


K)  cos0Rn  vj  dSQ  dS 


P 


(2.21) 


The  NxN  matrices  B  and  C  are  full  matrices  of  rank  N,  and  are  therefore  invertible 
(see  Section  2.6).  They  are  most  easily  constructed  if  v  corresponds  to  the  assumption  of  a 
constant  field  over  each  element  and  w  corresponds  to  collocation  at  centroidal  nodes  [see,  e.g., 
DeRuntz  and  Geers,  1978].  The  elements  of  B  and  C  are  then  given  by 
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(2.22) 


where  -  |  R;j  |  is  the  distance  from  the  centroid  of  element  i  to  an  integration  point  in 

element  j,  5-  is  the  Kronecker  delta,  and  is  the  angle  between  Rjj  and  the  surface  normal 
(going  into  the  fluid)  at  an  integration  point  in  element  j. 

To  obtain  the  semi-discretized  form  of  (2.17),  i.e..  that  produced  when  (2.17)  is 
discretized  in  space  but  not  in  time,  we  simply  premultiply  (2.20)  through  by  B“\  which 
yields 

DAA,  (t):  p(t)  +  cB"lCp(t)  -  pcu(t)  (2.23) 

It  would  be  advantageous  if  this  relation  were  converted  into  a  symmetric  form.  We 
accomplish  this  by  first  discretizing  (2.1 1)  to  obtain  the  matrix  LTA, 

B-‘Cp(t)  -  pu(t)  (2.24) 

Next,  we  obtain  a  suitable  boundary-element  expression  for  the  kinetic  energy  of  an 
inviscid,  incompressible  fluid  undergoing  irrotational  flow  (see  the  last  paragraph  in  Section 
2.3).  We  start  with  the  known  continuum  expression  (see.  e.g.,  Milne-'Piomson,  1960) 


T(t)  -  |  J  fcpfl)  uP(t)  dSP 


(2.25) 


where  the  asterisk  over  pP(t)  denotes  a  time  integration.  Introducing  the  discretization 
expressions  (2. 1 8).  we  then  obtain 


T(t)  -  J  iiT(t)A$(t) 


(2.26) 
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where  the  generalized  area  matrix  A  is  given  by 


A 


VQ  dSq 


(2.27) 


Note  that  A  is  a  diagonal  matrix  of  element  areas  if  v  corresponds  to  the  assumption  of  a 
constant  field  over  each  element. 

Now,  by  (2.24),  p(t)  and  u(t)  cannot  both  be  free  vectors;  hence  we  choose  u(t)  as  the 
free  vector  and  employ  (2.24)  to  introduce  p(t)  -  pC“'Bu  into  (2.26),  which  yields 

T(t)  -  J  puT(t)  AC'B  u(t)  (2.28) 

Then  we  separate  the  matrix  AC^B  into  its  symmetric  and  skew-symmetric  parts  and  note 
that  the  latter  contributes  nothing  to  the  kinetic  energy.  Hence  (2.28)  becomes 

T(t)  -  i  piiT(t)  (AC-'B)  ti(t)  (2.29) 

where  the  angular  brackets  denote  symmetrization  of  the  matrix  within. 

As  the  next  step,  we  treat  p(t)  as  a  prescribed  vector  and  write  the  fluid  work- 
potential  expression 


II(t)  «  pP(t)  uP(t)  dSP  -  uT(t)Ap(t)  (2.30) 

Js 

Thus,  on  the  basis  of  (2.29)  and  (2.30),  Hamilton’s  Principle,  5j(T+n)dt  »  0,  applied  for 
variations  Su,  yields 


p<AC"‘B)u(t)  -  Ap(t)  (2.31) 

The  matrix  p<AC“‘B)  is  known  as  the  fluid  mass  matrix  (DeRuntz  and  Geers,  1978),  which  is 
a  generalized  form  of  Lamb’s  inertia  coefficients  for  hydrodynamic  flow  about  rigid  bodies 
(Lamb,  1945).  The  fluid  mass  matrix  is  positive-definite. 


12 


For  our  purposes,  we  reverse  (2.31)  and  then  multiply  through  by 
A<AC“lB)_1  to  obtain 


A(AC-1B>-1Ap(t)  -  pAu(t) 


(2.32) 


But  because  A  is  symmetric, 

A<AC“1B)“1A  -  At<AC-,B>“1A 

-  <AT[AC*‘B]'1A> 

-  (AT^CA-'A)  (2.33) 

-  <AB->C> 

Hence  the  symmetric-matrix  LTA,  becomes  [cf.  (2.24)] 

<A  B-'C)  p(t)  -  pAu(t)  (2.34) 

A  way  to  get  this  result  more  directly  is  to  use  as  a  kinetic-energy  expression 
equivalent  to  (2.26) 


T(t)  -  |pT(t)Au(t)  (2.35) 

and  to  choose  p(t)  as  our  free  vector.  Employing  (2.24),  we  introduce  u  -  p-'B^Cp  into  (2.35) 
and  retain  only  the  symmetric  part  of  AB^C  to  obtain 

T(t)  -  i  p'‘pT (t)  (AB-‘C>  fct)  (2.36) 

Next,  we  treat  u(t)  as  a  prescribed  vector  and  write  for  the  fluid  work  potential 

n(t)  -  pP(t)  Up(t)  dS  -  pT(t)Au(t)  (2.37) 

Js 


Then  the  application  of  Hamilton's  Principle  for  variations  5p,  followed  by  double 
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differentiation  of  the  resulting  equation  in  time,  yields  (2.34). 

With  (2.34)  as  our  symmetric-matrix  LTA,.  our  symmetric-matrix  DA  A,  is  clearly  [cf. 

(2.23)] 

(DAA> ,(t):  Ap(t)  +  c<AB-'C)p(t)  -  pcAu(t)  (2.38) 

Using  (2.33),  we  can  write  this  equation  in  the  form 

Mp(t)  +  pcAp(t)  -  pcMu(t)  (2.39) 

where  M  -  p<AC~1B)  is  the  fluid  mass  matrix,  discussed  after  (2.31).  This  is  the  original 
form  of  DAA^  derived  by  inspection  in  Geers,  1971. 

2.6  MODAL  ANALYSIS  OF  THE  EXTERNAL  DAA,. 

Consider  the  following  eigenproblem  on  the  surface  S: 

c0"‘7  0q  -  X^P  (2.40) 

This  eigenproblem  pertains  to  Laplace’s  equation  for  the  conservative  problem  of  irrotational 
sloshing  of  an  inviscid,  incompressible  external  fluid  (see  the  last  paragraph  in  Section  2.3). 
Furthermore,  the  sloshing  problem  is  derivable  from  the  kinetic  energy  expression  for  the 
external  fluid,  which  is  a  positive  quadratic  form.  Hence  the  eigenvalues  Xn  are  real  and 
positive,  and  the  eigenfunctions  \JiPn  are  real  and  possess  the  orthogonality  property 

f  ^Pm^Pn  dSP  -  An5mn  (2.41) 

Js 

where  An  is  a  normalization  constant  and  5mn  is  the  Kronecker  delta. 

Following  standard  modal  analysis  procedures,  we  expand  pP(t)  and  up(t)  as 
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(2.42) 


oo 

Pp(t)  -  *p»  Pnw 

n-o 

00 

upW  -  21  ^Pn  u"(t) 
n-o 

introduce  them  into  (2.17),  employ  (2.40),  multiply  through  by  \im,  integrate  over  S.  and  utilize 
(2.41)  to  obtain  the  modal  DAAX  equations  for  an  external  acoustic  domain 

DAA"(t):  pft  +  XnPn  -  pcii,,  (2.43) 

This  result  shows  that  the  fluid  boundary  modes  for  Laplace's  equation  in  the  external  domain 
can  be  used  to  decompose  the  DAAt  into  uncoupled  modal  equations  (Geers,  1978). 

Modal  analysis  of  the  unsymmetric-matrix  DAAt.  (2.23).  proceeds  in  similar  fashion. 
The  pertinent  eigenproblem  is 


cB -‘C*  -  X# 


(2.44) 


and  the  orthogonality  statement  is 


"  5. 


mn 


(2.45) 


We  expand  p(t)  and  u(t)  as 


oo 

P(t)  -  YL  PnW 
n-o 

(2.46) 

00 

U(t)  U„(t) 

n-o 

introduce  them  into  (2.23),  employ  (2.44),  premultiply  through  by  and  utilize  (2.45)  to 
obtain  (2.43). 
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Modal  analysis  of  the  symmetric-matrix  DAAi,  (2.38),  differs  only  slightly  from  that  of 
its  unsymmetric  counterpart.  Instead  of  (2.44),  the  pertinent  eigenproblem  is 

c<AB-'C>*  -  XA*  (2.47) 


and  the  orthogonality  statement  is 


*£A*n-An«mn  (2.48) 

Proceeding  as  before,  we  introduce  the  modal  expansions  (2.46)  into  (2.38),  employ  (2.47), 
premultiply  through  by  and  utilize  (2.48)  to  obtain  (2.43). 

Although  the  continuum-operator,  unsymmetric-matrix.  and  symmetric-matrix  DAA/s 
all  produce  (2.43),  the  three  sets  of  modes  all  differ  slightly  from  one  another,  depending  upon 
the  choice  of  shape  functions  vP  and  weight  functions  wp.  and  the  degree  of  surface  mesh 
refinement.  The  numerical  determination  of  fluid  boundary  modes  for  surfaces  of  general 
geometry  is  discussed  by  DeRuntz  and  Geers.  1978. 
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SECTION  3 


DAA,  FOR  AN  INTERNAL  ACOUSTIC  DOMAIN 

Development  of  the  first-order  DAA  for  an  internal  acoustic  domain  is  complicated  by 
the  existence  of  low-frequency  dilatational  motion,  which  does  not  occur  in  an  external 
domain.  Hence,  while  ETA,  is  clearly  the  same  for  both  internal  and  external  domains,  LTA, 
and  thus  DAA,  for  the  internal  domain  differ  from  their  external  counterparts. 

3.1  EQUIVOLUM1NAL  AND  DILATATIONAL  FIELDS  AT  LOW  FREQUENCIES. 

We  recall  the  conservation-of-mass  equation,  the  constitutive  equation,  and  the  small- 
perturbation  assumption  for  an  acoustic  fluid  (see,  e.g..  Pierce.  1981) 

|f  +  V  (pu)-0  .  If-c’Jfi  .  V-(pu)2spV  u  (3.1) 

These  may  be  combined  to  yield 

V  u  -  -  (pc2)-1  |£  (3.2) 

In  order  to  accommodate  dilatational  fluid  motion  in  the  internal  domain,  we  take  p(r,t)  - 
pd(t),  so  that  (3.2)  becomes,  after  integration  in  time, 

V-u(r.t)  -  -  (pc2)-1  pd(t)  (3.3) 

Note  that  pd(t)  -  0  for  an  external  medium,  in  order  that  the  boundary  condition  of  zero 
acoustic  pressure  at  infinity  may  be  satisfied. 

Now  (3.3)  is  an  equation  that  holds  at  every  point  in  the  fluid  volume.  Hence  we  may 
integrate  it  over  the  volume  and  apply  the  divergence  theorem  to  the  left  side  of  the  resulting 
equation  to  obtain 

u(S,t)  dS  -  (pc2)-1  Vpd(t)  (3.4) 

Js 
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where  we  recall  from  Section  2.1  that  u  -  un  and  where  V  is  the  volume  of  the  internal 
fluid  domain.  Let  us  investigate  the  nature  of  the  solutions  to  (3.4). 

First,  we  take  the  fluid-particle  displacement  field  as  comprised  of  two  parts:  u(r.t)  - 
ue(r,t)  +  ud(r,t),  where  u®  (r.t)  is  the  homogeneous  solution,  i.e.,  that  for  which  pd(t)  -  0.  and 
ud(r,t)  is  the  particular  solution  produced  by  pd(t).  On  this  basis.  (3.4)  yields 

u®(S,t)  dS  -  0 
JS 

(3.5) 

ud(S,t)  dS  -  (pc2)-1  Vpd(t) 

Js 

We  recognize  ue  and  ud  as  linearly  independent  equivoluminal  and  dilatational  fluid-particle 
displacement  fields,  respectively.  Similarly,  p(r,t)  -  pe(r.t)  +  pd(t).  where  pe(S,t)  satisfies  a 
zero-average  equation  like  the  first  of  (3.5). 

Next,  we  show  that  ud  is  constant  over  the  surface  S  and  determine  its  relationship  to 
pd.  Suppose  that 


ud(S,t)  -  ud(t)  +  uv(S,t)  (3.6) 

in  which  ud(t)  is  the  average  of  ud(S,t),  given  by  ud(t)  -  aud(S.t).  where  a  is  the  averaging 
operator  defined  as 


«Qq  e 


(3.7) 


in  which  A  is  the  area  of  the  surface  S.  Then  integration  of  (3.6)  over  S  yields 


uv(S,t)  dS  -  0  (3.8) 

Js 
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But  if  uv  satisfies  this  equation,  then,  from  the  first  of  (3.5),  it  must  be  part  of  the 
homogeneous  solution  ue  rather  than  part  of  the  particular  solution  ud.  Hence  ud(S,t)  -  ud(t). 
and  the  second  of  (3.5)  yields 


pd(t)  -  pcT(A/V)ud(t)  (3.9) 

3.2  FIRST-ORDER  LATE-TIME  APPROXIMATION:  LTA,. 

Determining  LTA,  for  an  internal  acoustic  domain  requires  separate  consideration  of 
equivoluminal  motion  and  dilatational  motion.  For  equivoiuminal  motion,  where  the  flow  is 
incompressible,  LTA,  is  determined  in  the  same  manner  as  that  used  for  general  motion  in 
the  external  domain,  and  is  given  by  [cf.  (2.10)] 

/S~*7  Pq(s)  -  p  s2  Up(s)  (3.10) 

Note,  however,  that  7  here  pertains  to  an  inward  normal,  while  7  in  (2.10)  pertains  to  an 
outward  normal. 

The  preceding  equation  is  not  valid  for  dilatational  motion,  because  7Pq  -  0  when  pq 
is  constant  over  S.  This  is  readily  shown  by  introducing  into  Green's  second  identity 

j^P^q  -  qV*p)  dV.|s[p|a-q|2jdS  (3.1 1) 


the  particular  functions  p  -  1  and  q  -  1/RpQ,  the  latter  being  the  singular  solution  to 
Laplace's  equation.  In  this  case,  all  terms  vanish  except  that  produced  by  the  first  integrand 
on  the  right,  yielding 


{sA(l/RpQ)dS- 

-  -  |^R^J  cos^Rn  dS  -  -71  -0 


(3.12) 
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where  we  have  used  the  second  of  (2.9). 

To  determine  LTAj  for  dilatat ional  motion,  we  take  pq(s)  -  pd(s)  and  Uq(s)  -  ud(s). 
introduce  (2.7)  into  (2.4).  and  retain  on  both  sides  all  terms  through  those  of  order  s2  to  obtain 

|o  -  ^(s/c^tjl  j  pd(s)  -  ps*/S  1  ud(s)  (3.13) 

where  the  spatial  operator  jj  is  defined  as 


9Qq  3 


I 


eos^Rn  Qq  dSq 


(3.14) 


Because  (3.13)  must  agree  with  (3.9).  ijl  -  -2(V/A)01. 

To  derive  a  first-order  LTA  for  general  motion  in  an  internal  domain,  one  might 
simply  introduce  (2.7)  into  (2.4).  retain  on  both  sides  all  terms  through  those  of  order  s2,  and 
premultiply  through  by  fS~l  to  get 

j8"‘7PQ(s)  -  |(s/c)2/5-1tjpq(s)  -  ps2uP(s)  (3.15) 

For  equivoluminal  motion,  this  expression  contains  on  the  left  both  0(s°)  and  0(s2)  terms; 
hence  it  is  not  a  first-order  LTA.  The  correct  first-order  LTA  is 

LTA,(s):  P~11Pq(s)  +  (s/c)2LapQ(s)  -  ps2uP(s)  (3.16) 

Where  L  -  V/A.  For  equivoluminal  motion,  this  relation  becomes  (3.10)  because  ctpq(s)  -  0; 
for  dilatational  motion,  it  corresponds  to  (3.9)  because  7pq(s)  -  7Pd(s)  -  0.  apq(s)  -  a pd (s)  « 
pd(s),  and  Up(s)  -  ud(s).  We  note  that  (3.16)  is  spatially  non-local  and  singly  asymptotic  in  the 
limit  s  -+  0. 

3.3  FIRST-ORDER  DOUBLY  ASYMPTOTIC  APPROXIMATION:  DAA,. 

We  seek  here  an  expression  that  naturally  reduces  to  ETA,.  (2.5),  and  LTA,,  (3.16),  in 
the  limits  s  -*•  x>  and  s  -*  0,  respectively.  One  that  does  so  is 
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CsJ+3pP(s)  +  ^~17Pq(s)  +  (s/c^LapgCs)  -  CpcsJ+4uP(s)  +  ps2uP(s) 


(3.17) 


where  C  is  an  arbitrary  constant  and  j  ^  0.  As  discussed  in  Section  2.4,  however,  the 
indeterminacy  of  C  and  j,  and  the  presence  of  high-order  derivatives  prompts  us  to  reject  it 
as  a  first-order  DAA. 

Having  gained  in  Section  2.4  some  experience  in  the  method  of  operator  matching,  we 
propose  a  DAA,  trial  equation  that  immediately  satisfies  (2.5)  for  s  -*■  oo,  viz.. 

spP(s)  +  cPoPqts)  -  pcs*uP(s)  +  pc*s£/,Uq(s)  (3.18) 

where  P0  and  Ux  are  unknown  spatial  operators. 

For  equivoluminal  motion  with  s  •*  0,  we  write  (3.18)  as 

[Po  +  0(s)]pq(s)  -  pcs£/,Uq(s)  +  pslup(s)  (3.19) 

The  only  way  in  which  this  can  match  (3.16)  with  pQ(s)  -  Pq(s)  is  if  Ux  -  F,a.  which 
annihilates  the  first  term  on  the  right  side  of  (3.19).  and  if  P0  -  p~ly. 

Next,  we  consider  (3.18)  as  it  would  apply  to  dilatational  motion  with  s-K).  Dividing 
through  by  s,  and  noting  that  pP(s)  -  pd(s).  P0Pq($)  -  ^"17Pq(s)  -  0,  and  i/,Ug(s)  »  F,o!Uq(s)  - 
F,ud(s),  we  obtain 


pd(s)  -  pc2F,ud(s)  +  0(s)  (3.20) 

This  matches  (3.16)  with  pQ(s)  -  pd(s)  only  if  Vx  -  L_1. 

The  introduction  of  these  results  into  (3.18)  produces  the  first-order  DAA  for  an 
internal  acoustic  domain,  expressed  in  transform  space  as 

DAA|(s):  spP(s)  +  c^~‘7Pq(s)  -  pcs2uP(s)  +  pc^L^^auqls)  (3.21) 

and  in  the  time  domain  as 

DAA,(t):  Pp(t)  +  cj3-‘7pQ(t)  -  pcuP(t)  +  pc2L_1o(UQ(t)  (3.22) 

As  expected,  this  is  a  spatially  non-local  approximation. 
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3.4  MATRIX  DA  A,  FOR  BOUNDARY  ELEMENT  ANALYSIS. 


Petrov-Galerkin  semi-discretization  may  be  applied  to  the  internal  DA  A,  in  the  same 
manner  in  which  it  was  applied  in  Section  2.5  to  the  external  DAA,.  The  DA  A,  equation 
(3.21)  is  preoperated  through  by  /?.  the  discretization  formulas  (2.18)  are  inserted,  the 
weighted-residual  procedure  is  implemented,  and  the  resulting  matrix  equation  is 
premultiplied  through  by  B-1  to  yield 

DAA,(s):  sp(s)  +  cB"'Cp(s)  -  pcs2u(s)  +  pc2 V"‘s  1  aT  u(s)  (3.23) 

where  B  and  C  are  given  by  (2.21),  1  is  the  unit  vector,  and 


a  - 


1 


vQ  dSQ 


(3.24) 


Note  that,  because  VqI  «  1  [see  the  discussion  following  (2.18)].  aT  1  -  A.  Equation  (3.4.1) 
becomes  in  the  time  domain 

DAA, (t):  p(t)  +  cB"*Cp(t)  -  pcu(l)  +  pc2 V*1 1  aT  ii(t)  (3.25) 

Because  7pd(t)  -  0  and  7Pq(t)  *  0.  Cpd(t)  -  Clpd(t)  -  0  and  Cpe(t)  #  0.  Hence  the 
NxN  matrix  C  is  of  rank  N-l,  transforming  p(t)  into  a  vector  in  the  RNl  subspace  of 
equivoluminal  solutions.  Similarly,  because  aUg(t)  ■>  0  and  a  ud  (t)  #  0.  aT  ue(t)  -  0  and 
aTud(t)  «  aT  1  ud  (t)  -  Aud(t)  *  0;  hence  the  NxN  rank-one  matrix  laT  transforms  u(t)  into  a 
vector  in  the  Rl  subspace  of  dilatational  solutions.  As  in  the  case  of  external  fluid,  B  is  of 
rank  N. 

Here  too  B  and  C  are  most  easily  constructed  if  v  corresponds  to  the  assumption  of  a 
constant  field  over  each  element  and  w  corresponds  to  collocation  at  centroidal  nodes.  The 
elements  of  B  and  C  are  then  given  by  (2.22),  and  a  is  a  vector  of  element  areas. 

It  is  possible  to  symmetrize  (3.23)  and  (3.25).  To  do  so,  we  first  consider  (3.23)  for 
late-time  equivoluminal  motion  by  taking  p(s)  -  pe(s)  and  u(s)  -  ue  (s)  with  s  ■*  0;  the  inverse 
Laplace  transform  of  the  resulting  equation  is  [cf.  (2.24)] 
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B-'CPeft)  -  piU t) 


(3.26) 


Then,  using  the  arguments  given  after  (2.34),  we  conclude  that  Hamilton's  Principle  yields 
instead  of  (3.26)  [cf.  (2.34)] 


(AB-’Qp.W  -  pAu^t) 


(3.27) 


is 


Next,  we  note  that  the  potential  energy  associated  with  dilatational  motion  of  the  fluid 


U(t)- 


(3.28) 


which,  after  introduction  of  the  discretization  expressions  (2.13),  becomes 

U(t)-  iuJ(t)Apd(t)  (3.29) 

But,  for  late-time  dilatational  motion,  i.e.,  for  p(s)  -  pd(s)  and  u(s)  -  ud(s)  with  s  •*  0,  (3.23) 
yields 

p^sJ-peV-UaTu^s)  (3.30) 

Inverse  Laplace  transforming  this  equation  and  introducing  the  result  into  (3.29),  we  obtain 

U(t)  -  ^  pcJV_1ud(t)(A  laT>ud(t)  (3.31) 

where  only  the  symmetric  part  of  AlaT  has  been  retained. 

Next,  we  treat  pd(t)  as  a  prescribed  vector  and  write  the  fluid  work-potential 
expression 

n(t)  -  p£(t)  Up(t)  dSP  -  ud(t)Apd(t)  (3.32) 

Js 
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Then,  using  (3.31)  and  (3.32).  we  apply  Hamilton's  Principle  to  obtain 

Apd(t)  -  pc2  V“l(A  laT>ud(t)  (3.33) 

From  (3.25),  (3.27).  and  (3.33).  we  conclude  that  the  symmetric  counterpart  to  (3.4.3)  is 

<DAA),(t):  Ap(t)  +  c(AB*‘C>p(t)  -  pcAu(t)  +  pc*V-,<A  1  aT>ii(t)  (3.34) 

For  the  spatial  discretization  scheme  in  which  v  corresponds  to  the  assumption  of  a  constant 
field  over  each  element  and  w  corresponds  to  collocation  at  centroidal  nodes,  the  elements  of 
the  diagonal  matrix  A  and  of  the  vector  a  are  merely  the  areas  of  the  finite  elements,  and 
AlaT  is  already  symmetric. 

3.5  MODAL  ANALYSIS  OF  THE  INTERNAL  DAA,. 

Consider  for  the  internal  domain  the  eigenproblem  given  over  the  surface  S  by 

c/r‘70Q  -  Xtfp  (3.35) 

One  solution  to  this  eigenproblem  pertains  to  dilatational  motion,  for  which  \Jiq  -  yird  (constant 
over  S),  7^d  -  0.  and  so  Xd  -  0.  The  remaining  eigenvalues  and  eigenfunctions  pertain  to 
Laplace's  equation  for  the  conservative  problem  of  irrotational  sloshing  of  an  inviscid. 
incompressible  internal  fluid.  These  equivoluminal  modes  derive  from  the  kinetic  energy  of 
the  internal  fluid  and  hence  possess  real,  positive  eigenvalues  and  real  eigenfunctions,  as  well 
as  the  orthogonality  property  (2.41).  Orthogonality  extends  to  the  dilatational  mode  if  we 
assign,  say,  the  zero  index  to  ^  and  note  that  (2.41)  becomes  the  first  of  (3.5)  if  either  of  the 
modal  ^-subscripts  in  (2.41)  is  zero. 

In  summary,  then,  the  spectrum  of  the  eigenproblem  (3.35)  is  an  infinite  set  of  discrete 
eigenvalues,  the  first  of  which  is  zero  and  the  rest  of  which  are  positive.  The  eigenfunction 
corresponding  to  the  zero  eigenvalue  is  a  real  constant  over  the  surface  S,  which  constitutes 
the  dilatational  mode.  The  rest  of  the  eigenfunctions  are  real  with  zero  average  value,  thus 
constituting  the  equivoluminal  modes. 

From  the  preceding,  expansions  given  by  (2.42)  may  be  introduced  into  (3.22),  and  the 
orthogonalization  process  described  after  (2.42)  may  be  employed  to  obtain  the  modal  DAAX 
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equations  for  an  internal  acoustic  domain 


p0  -  pcu0  +  pc2!--1^ 

DAA|(t):  (3.36) 

Pn  +  Xn  Pn  -  PCK  .  n  ^  1 

These  equations  demonstrate  that  the  internal  DAA,  may  be  decomposed  into  uncoupled 
modal  equations  through  the  use  of  fluid  boundary  modes,  the  first  equation  pertaining  to  the 
dilatational  mode  and  the  second  to  the  equivoluminal  modes. 

Modal  analysis  of  the  unsymmetric-matrix  DAA,  is  similarly  straightforward.  The 
pertinent  eigenproblem  is  written  as  (2.44).  and  yields  eigenfunctions  Xn  and  eigenvectors  ^n. 
The  first  eigenvalue  Xe  «  0  and  the  corresponding  eigenvector  is  merely  the  unit  vector  1 
multiplied  by  a  normalization  constant;  the  remaining  N-l  eigenvectors  are  equivoluminal. 
The  orthogonality  statement  is  (2.45)  and  the  modal  expansion  is  given  by  (2.46).  Application 
to  (3.25)  of  the  orthogonalization  process  described  after  (2.46)  then  yields  (3.36),  inasmuch  as 
the  second  term  on  the  right  in  (3.25)  yields  *£laT>n  -  A5n0, 

Modal  analysis  of  the  symmetric-matrix  DAA,  follows  in  like  fashion.  We  know  that 
<AlaT)^n  vanishes  unless  m  -  n  -  0.  We  also  know  that  -  lc  1.  where  1„  is  a 
normalization  constant,  and  that  the  skew-symmetric  part  of  an  unsymmetric  matrix 
contributes  nothing  to  a  generalized  inner  product;  hence  ^J(AlaT>^0  -  ^ J A 1  aT  1D 1  - 
^jA0oaT  1  -  A0A.  where  we  have  used  (2.48)  and  the  statement  following  (3.24).  Therefore, 
application  to  (3.34)  of  the  orthogonalization  process  described  after  (2.48)  yields  (3.36). 
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SECTION  4 


DAA,  FOR  AN  EXTERNAL  ACOUSTIC  DOMAIN 

In  the  previous  sections  we  dealt  with  first-order  DAA’s  for  both  external  and  internal 
domains.  In  this  section,  we  derive  a  second-order  DAA  for  an  external  domain  by  the  method 
of  operator  matching. 

4.1  SECOND-ORDER  EARLY-TIME  APPROXIMATION:  ETA2. 

The  second-order  early-time  approximation  (Felippa,  1980)  is,  in  transform  space,  [cf. 

(2.5)] 

ETA2(s):  (s-K:K^)pp(s)=pcs2u1,(s)  (4.1) 

where  Kp  is  the  local  mean  curvature  of  the  surface.  Like  ETA„  ETA2  is  a  local  approximation 
in  space  that  has  each  element  of  the  surface  S  independently  generating  a  curved  wave  that 
propagates  outwardly  into  the  fluid.  Hence  (4.2)  is  often  referred  to  in  the  literature  as  the 
curved  wave  approximation. 

4.2  SECOND-ORDER  LATE-TIME  APPROXIMATION:  LTA2. 

The  derivation  of  LTA2  consists  merely  of  extending  that  of  LTA,,  i.e.,  introducing  (2.7) 
into  (2.4)  and  retaining  terms  of  order  s°  and  s1  on  the  left  and  s2  and  s3  on  the  right.  This  yields 

J  RkJ  cos^  pQ(s)  dSQ  =  pjs  2(Rp^  -  c  s)  Uq(s)  dSQ  (4.2) 

s  s 

Note  that  the  term  of  order  s1  on  the  left  has  vanished  identically. 

With  the  spatial  operator  definitions  (2.9)  and  (3.7),  we  can  express  (4.2)  in  operator  form 


LTA2(s): 


P_1Y  Pq(s)  =  ps 2  [up(s)  -  C'sA  p'‘a  uQ(s)] 


(4.3) 


4.3  SECOND-ORDER  DOUBLY  ASYMPTOTIC  APPROXIMATION:  DAA2. 

Pursuing  a  procedure  similar  to,  but  more  complicated  than,  that  for  DAA„  we  assume 
the  DAA2  trial  equation 

s^s)  +  csP ,  pQ(s)  +  c  2P0  pQ(s)  =  pcs3up(s)  +  pc  2s  2\J2  uQ(s)  (4.4) 

where  P„,  P,  and  U2  are  spatial  operators  to  be  defined.  For  s— >0,  we  write  this  as 

[Po  +  c  ’’s  Pj  +  0(s  2)]pQ(s)  =  ps  2(U2  +  c  ',s)uQ(s)  (4-5) 

where  we  adopt  the  convention  that  a  scalar  multiplying  Uq(s)  [or  pQ(s)J  yields  the  product  of 
that  scalar  and  uP(s)  [or  pP(s)].  For  s-*»,  we  divide  (4.4)  through  by  s2  and  write  the  result  as 

[1  +  cs  ‘‘P,  +  0(s'2)]pQ(s)  =  pc(s  +  cU2)uQ(s)  (4.6) 

In  order  to  match  (4.5)  to  (4.3)  and  (4.6)  to  (4.1),  we  need  to  invert  the  operator 
ensembles  either  on  the  left  or  right  of  (4.5)  and  (4.6).  This  leads  to  four  possible  solution 
procedures,  which,  as  might  be  expected,  are  equivalent.  For  example,  let  us  invert  the  operator 
on  the  left  side  of  (4.5)  to  obtain,  for  s  — >  0, 

Pp(s)  =  ps2[l  -  c-'sP/’P,  ♦  CKs2)]?;1^  +  c-'s)uQ(s)  (4.7) 

inasmuch  as  [1  -  c^sP^P,  +  0(s 2)  ]  P0-1  [P0  +c''sP,  +  0(s2)]  =  1  +  0(s2).  Hence, 

keeping  terms  of  order  s°  and  s\  and  then  multiplying  through  by  the  operator  P  ’y,  we  obtain 
the  result 


P''Y  PQ(s)  =  ps  2p-'y  P;1  [U2  +  c  -'s(l  -  P.P/'U,)]  uQ(s)  (4-8) 

Matching  this  to  (4.3),  we  find 
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P'WU,  =  1 

p-'yP^P.p;1!!,  -  1)  «  Ap-‘a 


(4.9) 


Similarly,  we  invert  the  operator  on  the  right  side  of  (4.6)  to  obtain 

s  1  [  1  -  cs‘‘U2  +  0(s'2)][l  +  cs_1P,  +  0(s'2)]pQ(s)  =  pcup(s)  (4.10) 

inasmuch  as  s'*[l  -  cs'lU2  +  0(s‘2)](s  +  cU2)  =  1  +  0(s'2 ).  Hence,  retaining  terms  of  order  s'1  and 
s'2,  and  then  multiplying  through  by  s2,  we  obtain 

spp(s)  +  c(P,  -  U2)pQ(s)  =  pcs  2up(s)  (4.11) 

Matching  this  to  (4.1),  we  find 

P,  -  U2  =  Kp  (4.12) 

The  unknown  operators  P0,  Pt  and  U2  may  now  be  determined  by  solving  (4.9)  and  (4.12) 
simultanously.  This  yields 

p0  =  xP'1? 

P,  =  %  +  kp  (413) 

U2 

where 

X  -  (P_1Y  -  Kp)(l  -  AP^aP^y)"1  (4-14) 

The  introduction  of  these  results  into  (4.4)  produces  the  second-order  DAA  for  an  external 
acoustic  domain,  expressed  in  transform  space  as 

DAA^s):  s  2pp(s)  +  csKppp(s)  +  csxpQ(s)  +  c^P^yp^s) 

(4.15) 

■  pcs  2[sup(s)  +  cxuQ(s)] 

and  in  the  time  domain  as 
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DAAj(t): 


(4.16) 


PP(t)  +  ctcppp(t)  +  cxf>Q(t)  +  c^yp^t) 

=  pc[u^(t)  +  cxuQ(t)] 

This  result  agrees  with  that  of  Nicolas-Vullienne,  1989.  In  his  formulation  the  term 
AP  'aP  ’y  in  the  expression  for  x  [see  (4.14)]  is  replaced  by  AP'ly(l/27t)aP'tY.  This  is  due  to  a 
result  obtained  previously  by  Ohayon,  1983,  which,  in  the  present  context,  states  that  y'a  reduces 
to  (l/2;t)a.  Extending  an  idea  outlined  by  Felippa  in  1980,  Nicolas- Vullierme  was  the  first  to 
formulate  DAA’s  by  operator  matching,  making  use  of  the  Fourier  transform  and  the  method  of 
stationary  phase.  Although  his  formulation  works  well  for  non-concave  surfaces,  it  becomes 
cumbersome  for  concave  surfaces,  a  circumstance  that  is  disadvantageous  for  internal  domains. 
This  is  because  high-frequency  behavior  does  not  correspond  to  early-time  behavior  when  a  ray 
departing  from  one  surface  element  along  the  latter's  normal  intersects  another  surface  element. 
In  the  time  domain,  causality  prevents  such  early- time  interaction. 


4.4  MATRIX  DAA2  FOR  BOUNDARY  ELEMENT  ANALYSIS. 


We  now  apply  Petrov-Galerkin  semi-discretization  to  the  external  DAA2.  To  avoid  the 
discretization  of  inverse  operators  in  (4.15),  we  spatially  discretize  ETA2,  given  by  (4.1),  and 
LTA2,  given  by  (4.3),  and  then  use  the  method  of  matrix  matching  to  generate  the  matrix  DAA2. 

We  start  with  (4.1).  Introducing  the  Laplace  transforms  of  (2.18)  and  implementing  the 
weighted-residual  procedure  described  in  section  2.5,  we  obtain  the  matrix  ETA2  in  transform 
space 


ETAj(s):  (sj  +  cK)  p(r)  =  pcs  2J  u(s) 

where 

J  =  J  wpVpTd5p 

s 

K  =  JwpKpVp JdSp 

s 


(4.17) 


(4.18) 
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Applying  the  same  procedure  to  (4.3)  after  it  has  been  preoperated  through  by  j$,  we  find  the 
matrix  LTA2 

LTA20t):  C  p(s)  ■  p$2(B  -  c  "'s a  aT)u(s)  (4.19) 

where  B  and  C  are  given  by  (2.21),  a  is  given  by  (3.24),  and 

a  =  J  wpdSp  (4.20) 

s 


For  the  simplest  spatial  discretization  scheme,  i.e.,  that  for  a  constant  field  asumed  over 
each  element  and  collocation  at  centroidal  nodes,  a  becomes  the  unit  vector,  a  becomes  a  vector 
of  element  areas,  J  becomes  the  identity  matrix,  and  K  becomes  a  diagonal  matrix  of  local  mean 
curvatures. 

Following  the  procedure  used  in  the  method  of  operator  matching,  we  propose  here  the 
trial  equation  [cf.  (4.4)] 

(s2J  +  csP1  +  c2P0)p(s)  «  pcs2(sj  +  cU2)u(i)  (4.21) 

where  P„  P,  and  U2  are  unknown  matrices.  For  s-»0,  let  us  premultiply  this  equation  through 
by  C  times  the  inverse  of  (s2J  +  csP,  +  c2P0)  and  then  write 

(s2J  +  csP.  +  c2P J"1  =  [c2P  (I  +  c  's P0'P.  +  c-'s2?;^)]-1 

(4.22) 

=  c  -2[I  -  c  -'s  P^P,  +  0(s2)]P0'* 

to  obtain 


Cp (s)  =  ps2CP;,[UJ  +  c  -‘s(J  -  P.Po-'U,)  +  0(s2)]u(s) 
We  then  match  this  to  (4.20)  through  order  s',  which  gives 

cp;‘u2  =  b 

CP;‘(J  -  PjPo^Uj)  =  -a  a  T 


(4.23) 


(4.24) 
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For  s — >°°,  let  us  premutilply  (4.21)  through  by  J  times  the  inverse  of  (sj  +  cU2)  and  then 

write 


(sj  +  cUj)'1  =  [sj(l  +  cs-,J-,UJ-'  =  s'[l  -  cs-’J^U,  +  0(s2)]J'*  (4.25) 

to  obtain 

[sj  +  c(P,  -  Uj)  +  0(5_1)]p(i)  =  pcs 2  J  u(s)  (4.26) 

We  match  this  to  (4.17)  through  order  s°,  which  yields 

P,  -  Uj  =  K  (4.27) 

The  unknown  matrices  P0,  P,  and  U2  may  now  be  determined  by  solving  (4.24)  and  (4.27) 
simultanously  to  obtain  [cf.  (4.13)] 

Pc  =  JXB‘C 

P,  =  J  X  +  K  (4.28) 

Ua  =  JX 

where  [cf.  (4.14)] 

X  =  (B  *‘C  -  J  -‘K)  (I  -  B  -xa  a  TB  -’C)'1  (4*29> 

Substituting  (4.28)  into  (4.22),  we  obtain  the  matrix  DAA2  for  an  external  acoustic  domain, 
written  in  transform  space  as  [cf.  (4.15)] 


DAA2(s): 


[s2I  +  cs(X  +  J-‘K)  +c2XB'1C]p(s) 

=pcs2(sl  +  cX)u(s) 


(4.30) 


and  in  the  time  domain  as  [cf.  (4.16)] 


DAAj(t): 


(4.31) 


p(t)  +  c(X  +  J-‘K)p(t)  +  c2XB-‘Cp(t) 

=pc[u(t)  +  cX  u(t)3 

This  result  differs  somewhat  from  that  of  Felippa,  1980,  who  applied  his  matching  procedure  to 
a  scalar  representation  of  the  true  integral  formulation.  The  extension  to  a  matrix  form  was  done 
inferentially,  which  incorrectly  ordered  some  matrix  multiplications.  Also,  in  contrast  to  DAAj, 
DAA2  does  not  lend  itself  to  symmetrization. 

4.5  OTHER  FORMULATIONS. 

Matching  is  not  the  only  way  to  derive  higher-order  DAA’s.  In  1978,  Geers  formulated 
a  matrix  DAA2  for  an  external  acoustic  domain  on  the  basis  of  fluid  boundary  modes.  This 
procedure  introduced  a  free  parameter  in  each  modal  DAA2  equation  that  could  be  used  to 
optimize  its  accuracy  across  the  intermediate  frequency  region.  This  has  not  been  deemed 
necessary  for  external-domain  problems,  but  may  be  found  advantageous  for  internal  domains 
(Geers,  1990). 

In  a  modal  DAA  formulation,  the  fluid  boundary  modes  remain  uncoupled  across  the 
entire  frequency  range.  This  does  not  reflect  the  situation  when  an  exact  steady-state  acoustic- 
radiation  formulation  is  decomposed  by  fluid-boundary-mode  analysis;  in  this  circumstance,  the 
modes  are  coupled  across  the  intermediate  frequency  region.  An  examination  of  (4.15)  from  the 
viewpoint  of  (2.40)  reveals  that  matched  DAA2’s  do  not  submit  to  fluid-boundary-mode 
decomposition,  in  consonance  with  the  situation  characterizing  an  exact  steady-state  analysis. 

In  addition  to  matched  DAA  formulations  and  modal  DAA  formulations,  there  are 
symmetric  DAA  formulations  (Geers  and  Zhang,  1988).  A  symmetric  formulation  was  found  to 
be  essential  in  achieving  satisfactory  performance  by  a  first-order  DAA  in  electromagnetic 
scattering  problems.  In  the  context  of  acoustic  scattering,  a  symmetric  DAA  involves  augmenting 
(2.4)  with  an  equation  obtained  by  taking  the  partial  derivative  of  (2.4)  with  respect  to  the  surface 
normal.  Interesting  work  remains  to  be  done  here. 
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SECTION  5 


DAA2  FOR  AN  INTERNAL  ACOUSTIC  DOMAIN 

From  Section  3,  we  recall  that  the  derivation  of  DAA’s  for  an  internal  acoustic  domain 
is  complicated  by  the  existence  of  low-frequency  dilatational  motion.  Hence,  LTA2  and  DAA2 
for  the  internal  domain  differ  from  their  external  counterparts,  while  ETA2  is  the  same  for  both 
internal  and  external  domains. 

5.1  SECOND-ORDER  LATE-TIME  APPROXIMATION:  LTA2. 

LTA2  for  equivoluminal  motion  in  an  internal  domain  is  here  determined  by  the  procedure 
previously  utilized  for  general  motion  in  the  external  domain;  it  is  given  by  [cf.  (4.3)] 

P"*Y  Pq(s)  ■  ps2  uPe(s)  (51) 

where  y  pertains  to  an  inward  normal  and  the  inverse  of  y  does  not  exist.  Note  that  the  term 
corresponding  to  the  last  term  in  (4.3)  is  absent  here  because  auQ'(s)  *  0. 

To  determine  LTA2  for  dilatational  motion,  we  replace  pQ(s)  by  pd(s)  and  Uq(s)  by  ud(s) 
in  (2.4),  introduce  (2.7)  into  (2.4),  and  retain  on  both  sides  all  terms  through  those  of  order  s3. 
We  then  introduce  operator  symbols  to  represent  the  spatial  integrals  and  recall  that 
y  pd(s)  =  0  and  a  Uq(s)  =  ud(s)  to  obtain 

[  0  -  ifs/c)^  1  +  I(s/c)3C  1 1  p  d(s)  "  Ps  2  1  ~  (s/c)A]u  d(s)  (5-2) 

«  J 

where  the  spatial  operators  P  and  tj  have  been  defined  in  (2.9)  and  (3.14),  and  q  is  defined  as 

C  qQ  s  /  Rpq  cos^  qQ  dSQ  (5.3) 

s 

As  discussed  in  Section  3.2,  a  naive  derivation  of  LTA2  for  general  motion  in  an  internal 
domain  consists  merely  of  the  introduction  of  (2.7)  into  (2.4)  and  the  retention  on  both  sides  of 
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all  terms  through  those  of  order  s3.  This  yields, 

YPq(s)  -  i-Cs/c)^  pQ(s)  +  i.(s/c)3C  pQ(s) 

l  j  (5.4) 

=  ps2[puQ(s)  -  (s/c)AauQ(s)] 

Comparing  this  with  (5.1)  and  (5.2),  we  see  that  the  0(s2)  and  0(s3)  terms  on  the  left  side  are 
not  appropriate  for  equivoluminal  motion  but  are  needed  for  dilatational  motion;  hence  we 
introduce  the  averaging  operator  a  into  these  terms  to  obtain 


[Y  +  (s/c^pa  +  -I(s/c)3£alpQ(s) 


LTAj(s): 


3 . QV 

=  ps2[p  -  (s/c)Aa]u0(s) 


(5.5) 


where  we  have  used  the  fact  that  i'll  =  -2Lpl,  as  established  in  Section  3.2. 


5.2  SECOND-ORDER  DOUBLY  ASYMPTOTIC  APPROXIMATION:  DAA2. 

To  construct  a  trial  equation  for  internal  DAA2,  we  increase  the  order  of  the  trial  equation 
for  the  internal  DAAl5  (3.18).  Thus  we  try 

(s2  +  csP,  +  c2P0)pQ(s)  =  pcs(s2  +  csU2  +  c2U,)uq(s)  (5.6) 

where  P0,  P„  U,  and  U2  are  spatial  operators  to  be  found.  Determining  these  operators  requires 
separate  consideration  of  equivoluminal  motion  and  dilatational  motion,  accompanied  by 
extensive  matching. 

Let  us  first  examine  the  trial  equation  (5.6)  for  equivoluminal  motion.  For  s-»0,  the  only 
way  that  this  equation  with  pQ(s)  =  pQ*(s)  and  uQ(s)  =  uQ*(s)  can  match  (5.1)  is  if  U,=V,a;  hence 
the  trial  equation  (5.6)  for  s-»0  becomes  for  equivoluminal  motion, 

[P0  +  c"*sP,  +  0(s2)]pQe(s)  =  ps2(U2  +  c  _1s)  Uq(s)  (5-7) 

For  s — >°°,  (5.6)  becomes,  for  equivoluminal  motion. 
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[1  +  cs-'P,  +  0(s  '2)]  pQe(s)  =  pc(s  +  cU2)  Uq(s)  (5-8) 

Following  the  same  procedure  as  that  used  for  general  motion  in  an  external  domain  (see 
Section  4.3),  we  determine  the  spatial  operators  P0,  Pi  and  U2  by  matching  (5.7)  to  (5.1)  and  (5.8) 
to  (4.1);  the  results  are  [cf.  (4.13)] 

P„  -  (P-’Y  -  Kp^Y 

P1  -  p-»y  (5.9) 

U2  =  p-'Y  -  Kp 


Let  us  remind  ourselves  that  y  here  pertains  to  an  inward  normal  and  the  inverse  of  y  does  not 
exist. 

We  now  observe  that  the  preceding  development  for  equivoluminal  motion  would  be 
unaffected  if  Pc,  P,  and  U2  were  replaced  by  P0  +  QoOt,  Pi  +  Qi<x  and  U2  +  V2a,  respectively. 
Hence,  remembering  that  U,  *  Via,  we  write  the  updated  trial  equation 

[s2  +  cs(Pj  +  Q,a)  +  c2(P0  +  Qoa)  ]  pQ(s) 

*  pcs  [ s 2  +  cs(U2  +  V2a)  +  c 2 V,a]  uQ(s) 

where  P0,  Pi  and  U2  are  known,  and  Q„,  Qt,  Vj  and  V2  are  unknown. 

Let  us  now  examine  (5.10)  for  dilatational  motion.  With  pQ(s)  =  pd(s),  apd(s)  =  pd(s), 
Uq(s)  =  ud(s),  aud(s)  =  ud(s)  and  yl  =  0,  (5.10)  may  be  expressed  for  dilatational  motion  as  [note 
(5.9)] 


[s2  +  csQj  +  c  2QC  ]  p  d(s) 

=  pcs[s2  +  cs(V2  -  CKp )  +  c2V,]ud(s) 


(5.11) 


Unique  solutions  for  Q,,,  Vi  and  V2  can  be  found  in  the  R1  subspace  of  dilatational  solutions. 
Hence,  we  project  (5.11)  into  this  R1  subspace  by  preoperating  through  by  a,  which  leads  to  the 
conclusion  that  Q,,  Q„  V,  and  V2  must  be  scalars.  This  yields  the  scalar  equation 
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(s2  +  csQj  +  c2Q0) p  d(s)  *  pcs [s2  +  cs  (V2  -  cic)  +  c2V,)]u  d(s)  (5-12) 

where  k  is  the  average  curvature  over  the  surface  S.  For  future  use,  we  write  this  equation  for 
small  s  as 


(c2Qo  +  csQ,  +  s 2 ) p d(s)  =  pc2s[cV,  +  s(V2  -  ck)  +  0(s2)]ud(s)  (5-13) 

and  for  large  s  as 

[1  +  cs_1Q,  +  0(s'2)]pd(s)  =  pcs[l  +  cs~l(V2  -  cic)  +  0(s'2)]ud(s)  (5.14) 

In  order  to  determine  Qo,  Q„  V,  and  V2  by  second-order  matching,  we  must  also  project 
(5.5)  and  (4.1)  onto  the  R1  subspace  of  dilatational  solutions.  Hence,  with  pQ(s)  =  pd(s)  and  Uq(s) 
=  ud(s),  we  preoperate  (5.5)  through  by  a  to  obtain  LTA2  for  dilatational  motion 

[Lb  +  -I(s/c)z]p  d(s)  =  pc2[b  -  (s/c)A]ud(s)  (5.15) 

where  b  =  a|3l  and  z  =  a£l.  Then  we  do  the  same  to  (4.1)  to  obtain  ETA2  for  dilatational 
motion 

(1  +  cs"*K)p  d(s)  *  pcsud(s)  (5.16) 

We  now  perform  second-order  small-s  maching  of  (5.13)  and  (5.15).  First,  we  observe 
that  (5.13)  can  match  (5.15)  as  s— >0  only  if  Q„  =  0.  Then  we  multiply  (5.13)  through  by 
(cV,)'1  [1  -  (s/c)V,‘1(V2  -  k)]  to  obtain 


{vr'Q,  +  (s/c)vr'[i  -  v.-'q/v,  -  k)]  +  o(S2)}Pd(S) 

*  pc2[l  +  0(s2)]u  d(s) 


(5.17) 


Next,  we  multiply  (5.15)  through  by  b'l[l  +  (s/c)b  ‘A]  to  get 
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(5.18) 


[L  +  (s/c)b-*(V  +  -jz)  +  0(s2)]pd(s)  =  pc2[l  +  0(s2)]  u  d(s) 


Matching  these  two  equations  through  O(s),  we  find 

vr‘Qt  =  l 

v/ll  -  Vr'QjCVj-K)]  =  b-‘(V  +  JZ) 


(5.19) 


We  then  perform  second-order  large-s  matching  of  (5.14)  and  (5.16).  Multiplying  (5.14)  through 
by  1  -  cs‘‘(V2  -  k),  we  get 


[1  +  cs“‘(Q,  -  V2  +  ic)  +  0(s  *2)]  p  d(s)  =  pcs[l  +  0(s  *2)]  u  d(s)  (5.20) 
Matching  this  to  (5.16),  we  find 

Q,  -  V2  +  ic  =  ic  (5.21) 


Finally,  we  solve  (5.19)  and  (5.21)  simutaneously  to  obtain 

Q,  -  Ld 
V,  -d 

V2  =  Ld 


where 


d  = 


L2 


1  +  kL 
+  (V  +  j  z)/b 


(5.22) 


(5.23) 


Thus,  by  introducing  (5.9),  (5.22),  and  Q,  =  0  into  (5.10),  we  obtain  the  second-order  DAA  for 
an  internal  acoustic  domain,  expressed  in  transfom  space  as 
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DAA^s): 


(5.24) 


s^s)  +  cs(P_1y  +  Lda)pQ(s)  +  c2((5_1y  -  k^IJ-'yPqCs) 

=  pcs[(s2  -  cskp)up(s)  +  cs(3_1y  +  Lda)uQ(s)  +  c2dauQ(s)] 

and  in  the  time  domain  as 


DAAj(t): 


pp(t)  +  c( p_1Y  +  Lda) pQ(t)  +  c2(p',Y  -  k^P^YPqW 

=  pc[u^(t)  -  cKpiip(t)  +  c(p_1Y  +  Lda)iiQ(t)  +  c2dauQ(t)] 


(5.25) 


5.3  MATRIX  DAA2  FOR  BOUNDARY  ELEMENT  ANALYSIS. 

We  will  now  derive  a  matrix  form  of  the  second-order  DAA  for  internal  acoustic  domains 
by  Petrov-Galerkin  semi-discretization.  The  matrix  ETA2  is  given  by  (4.17).  The  matrix  LTA2 
is  obtained  by  introducing  the  Laplace  transforms  of  (2.18)  into  (5.5)  and  then  implementing  the 
weighted  residual  procedure  described  in  Section  2.5.  The  result  is 

[C  +  (s/c^LA-'BIa7  +  i(s/c)3  A^Zla^pfr) 

LTAj(s):  3  (5.26) 

*  p  s  2  [B  -  (s/c)a  a  T]  u(s) 

where  B,  C,  a  and  a  have  been  given  in  (2.21),  (3.24)  and  (4.21),  respectively,  and  where 

Z  =  //wpRpq  C0S<^Rn  vqT  dSQ  dSr  (5.27) 

s  s 

For  equivoluminal  motion,  aTpe  =  0  and  aTue  =  0,  so  (5.26)  yields  as  LTA2  for 
equivoluminal  motion 

Cpe(s)  =  ps 2  B  ut(s)  (5.28) 

For  dilatational  motion,  pd(s)  =  lpd(s)  and  ud(s)  =  lud(s);  also,  because  yl  =  0,  Cl  =  0.  Thus, 


with  aTl  =  A,  (5.26)  becomes 


[LB  1  +  I  (s/c)  Z  1]  pd(s)  =  pc 2  [B  1  -  (s/c)  A  a]  ud(s)  (5.29) 


At  this  point,  we  introduce  the  DAA2  trial  equation  [cf.  (5.6)] 

(s2J  +  csP,  +  c2Pc)p(s)  =  pcs(s2J  +  csU2  +  c2U,)u(s)  (5.30) 

where  P0,  P,  U,  and  U2  are  unknown  matrices.  To  match,  as  s— >0,  LTA2  for  equivoluminal 
motion,  i.e.,  (5.28),  U,  must  be  of  the  form  vaT  where  v  is  an  unknown  vector,  so  that  the  last 
term  on  the  right  in  (5.30)  will  vanish.  Then  the  procedure  used  in  Section  4.4  to  obtain  matrix 
DAA2  for  an  external  domain  may  be  applied  here  to  obtain  [cf.  (4.28)  and  (5.9)] 

P„  =  (J  B  ‘C  -  K)B  -*C 

P,  =  J  B‘C  (531) 

U2  =  J  B  *C  -  K 


Now  the  preceding  equivoluminal  development  would  be  unaffected  if  P0,  P,.  U,  and  U2 
were  replaced  by  P„  +  Q„laT,  P,  +  Q,laT,  V,laT  and  U2  +  V2laT,  respectively.  Thus  we  write 
the  updated  trial  equation 

[  s 2  J  +  cs(P  +  Q,laT)  +  c2(P0  +  Q„laT)]p(s) 

(5.32) 

*  pcs[s2J  +  cs(U2  +  V2laT)  +  c2V,laT]u(s) 


where  P0,  P,  and  U2  are  known  matrices,  and  Qo,  Qi,  V,  and  V2  are  unknown  scalars. 

Let  us  now  consider  dilatational  motion.  With  p(s)  =  lpd(s),  u(s)  =  lud(s),  aTl  =  A,  and 
Cl  =  0,  (5.32)  yields  [note  (5.31)  and  cf.  (5.11)] 


(s2Jl  +  csAQ,l  +  c2AQ0l)pd(s) 

=  pcs[s2Jl  +  cs(AV21  -  Kl)  +  c2AV,l]ud(s) 


(5.33) 


Premultiplication  through  by  N'llT,  where  N  is  the  size  of  the  discrete  system  (i.e.,  the  number 
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of  degrees  of  freedom),  then  yields  [cf.  (5.12)] 


(s2J  +  csAQ,  +  c2AQo)  pd(s) 

(5.34) 

=  pcs[s2J  +  cs(AV2  -  K)  +  c2AV,]ud(s) 

where  J  =  N  '1TJ1  and  K  =  N'11TK1.  For  future  use,  we  write  this  equation  for  small  s  as  [cf. 
(5.13)] 

(c2Qo  +  csQj  +  s2A_,J)p(s) 

(5.35) 

=  pc2s[cV,  +  s(V2  -  A_,K)  +  0(s2)]ud(s) 
and  for  large  s  as  [cf.  (5.14)] 

[J  +  cs  AQ,  +  0(s'2)]pd(s) 

(5.36) 

=  pcs[J  +  cs"‘(AV2  -  K)  +  0(s'2]ud(s) 

In  order  to  determine  our  unknown  dilatational  coefficients  by  matching,  we  must  also 
project  (5.29)  and  (4.17)  onto  the  R1  space  of  dilatational  solutions.  Hence,  with  p(s)  =  lpd(s) 
and  u(s)  =  lud(s),  we  premultiply  (5.29)  through  by  N  llT  to  obtain  LTA2  for  dilatational  motion 
[cf  (5.15)] 

[LB  -  i.(5/c)Z]^(5)  =  pc2[B  -  ( s/c)L-'Y]ud{s )  (5.37) 

where  B  =  N  lITBI,  Z  =  N  '1tZ1,  and  Y  =  N''VlTa.  Then  we  do  the  same  to  (4.4.1),  which 
yields  ETA2  for  dilatational  motion  [cf.  (5.16)] 

(J  +  cs_1K)pd(s)  =  pcsJud(s)  (5.38) 

We  are  now  ready  for  matching.  First,  we  note  that  (5.35)  can  match  (5.37)  for  s— >0  only 
if  Q„  =  0.  Then,  we  follow  the  same  procedure  as  that  employed  in  the  previous  section  to 
obtain  [cf.  (5.22)] 
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where  [cf.  (5.23)] 


Q,  =ld 


V,  -D 


V2  =  LD 


(5.39) 


(J  ♦  KL)/A 
L2  +  (Y  +  yZ)/B 


(5.40) 


Thus,  by  introducing  (5.31),  (5.39)  and  Q„  =  0  into  (5.32),  we  obtain  the  second-order  matrix 
DAA  for  an  internal  acoustic  domain,  expressed  in  transform  space  as  [cf.  (5.24)] 


DAAj(s): 


[s2I  +  cs(B‘C  +  LDJ-‘laT)  +  c2(B‘C  -  J  _,K)B  *‘C]  p(s) 

=  pcsjs 2I  +  cs(B-‘C  -  J-‘K)  +  LDJ *,laT]  +  c2DJ -‘la7 Ju(s) 


(5.41) 


and  in  the  time  domain  as  [cf.  (5.25)] 


DAA2(t): 


p(t)  +  c(B  "‘C  +  LDJ-'laT)^(s)  +c2(B-‘C  -  J  -*K)B  "‘C  p(t) 

=  pc [ u(t)  +  c(B-‘C  -  J-‘K  +  LDJ-‘laT)u(t)  +  c2DJ-‘laTu(t)] 


(5.42) 


In  contrast  to  DAA„  DAA2  does  not  lend  itself  to  symmetrization. 
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SECTION  6 


MODAL  EQUATIONS  FOR  A  SPHERICAL  GEOMETRY 

In  this  section,  modal  exact  and  DAA  equations  are  derived  for  axisymmetric  flow  in  a 
spherical  geometry,  in  which  the  wave  equation  separates,  thereby  admitting  solution  by 
separation  of  variables  (Morse  and  Ingard,  1968).  The  derivation  utilizes  velocity-potential  fields 
external  and  internal  to  a  spherical  surface  and  introduces  residual  potentials  to  facilitate  the 
development  (Geers,  1969,  1971,  1972).  Exact  modal  equations  linking  the  velocity  potentials 
and  their  derivatives  are  first  obtained,  from  which  early-time  and  late  time  approximations  are 
generated.  Then  modal  DAA’s  of  first  and  second  order  are  constructed  by  scalar  matching,  the 
scalar  form  of  operator  and  matrix  matching.  A  nondimensional  formulation  is  used  throughout 
the  section,  with  length  normalized  to  a,  the  radius  of  the  sphere,  time  normalized  to  a/c,  and 
pressure  to  pc2. 

6.1  EXACT  MODAL  EQUATIONS  FOR  THE  EXTERNAL  FLUID. 

For  axisymmetric  flow,  the  wave  equation  in  spherical  coordinates  admits  solutions  of  the 
form  (Morse  and  Ingard,  1968) 

M 

4>(r,e,t)  =  5>B(r,t)Pn(cos0)  (61> 

n*0 

where  r  and  0  are  the  radial  and  meridional  coordinates,  respectively,  <J>(r,0,t)  is  the  velocity 
potential,  and  P„(cos0)  is  the  Legendre  polynomial  of  order  n  (Abramowitz  and  Stegun,  1964). 
By  taking  the  Laplace-transform  of  both  sides  of  (6.1)  and  utilizing  the  orthogonality  property 
of  Legendre  polynomials,  one  can  obtain  the  following  ordinary  differential  equation  for  each  of 
the  <J)n(r,s): 

^-21  +  2 S-i  -  a2  -  n(n  +  l)K>n  =  0  (6.2) 

d£2  dq 

where  £  -  rs.  The  solution  to  this  equation  that  vanishes  as  ^  is 
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<t>„(r,s)  =  fn(s)  kn(rs) 


(6.3) 


where  fn(s)  is  an  underdetermined  function  and  k„(^)  is  the  nth-order  modified  spherical  Bessel 
function  of  the  third  kind  (Abramowitz  and  Stegun,  1964). 

Geers  (1969,1971,1972)  has  shown  that  <J>n  at  r  =  1  is  conveniently  obtained  as  the  solution 
of  the  equation 

<J>^(s)  +  s<j>n(s)  +  <J>n(s)  +  v^b(s)  =  0  (6.4) 

where  the  r-subscript  denotes  radial  differentiation,  underlining  indicates  location  on  the  surface 
r  =  1,  and  the  modal  residual  potential  \j/n(s)  is  given  by 

HL  (s)  -  -  [s  +  1  +  s^l]i(s)  (6-5) 

"  kn(s)  ■ 

in  which  k„’  is  the  derivative  of  k„  with  respect  to  its  argument  But  k„(£)  is  given  by 

(6.6) 

m-0 

where  rnm  =  (n+m)!  [2mm!  (n-m)!]'1.  Hence  (6.5)  and  (6.6)  yield 


Ern  nS“^(s)  =  EmrnraSn"niD(s)  (6-7) 

m^)  0i»O 

Equations  (6.4)  and  (6.7)  constitute  two  equations  for  the  two  unknowns  <j>„(s)  and  \j/„(s) 
in  terms  of  the  radial  derivative  4i^(s)-  They  are  key  equations  in  the  exact  fluid-structure 
interaction  formulation  of  Section  7. 


6.2  EXACT  MODAL  EQUATIONS  FOR  THE  INTERNAL  FLUID. 

We  now  apply  the  approch  of  the  preceding  section  to  the  internal  acoustic  fluid. 
Equation  (6.1)  and  (6.2)  carry  over,  but  (6.3)  become 
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<5>„(^)  =  ffs)  ifrs) 


(6.8) 


where  fn(s)  is  again  an  unknown  function  and  i„(<;)  is  the  nth-order  modified  spherical  Bessel 
function  of  the  first  kind,  which  remains  finite  at  r  =  0  (Abramowitz  and  Stegun,  1964). 

Preceding  as  before,  we  find  that  <J>„  at  r  =  1  is  conveniently  obtained  as  the  solution  of 
equation  [cf.  (6.4)] 

<j)^(s)  -  s<j>n(s)  +  <j>n(s)  +  )£a(s)  =  0  (6.9) 

where  the  modal  residual  potential  \&(s)  is  given  by  [cf.  (6.5)] 

i  *(s) 

V  (s)  -  [s  -  1  -  s__.]<|>  (s)  (6-10) 

■  i„(s)  " 

where  i„’  is  the  derivative  of  i0  with  respect  to  its  argument.  But  i„(£)  is  given  by 

in(£)  =  rnm4n'm ((-Dm  -  (-l)ne^]e^ 

2  m*0 

so  that  (6.10)  is  equivalent  to  the  delayed-differential  equation 

E(-1)"rnms"X(s)  =  ^(-irms—r^s) 

m*0  mM) 

n 

+  (~1)nIT  (»)  ”  2si  (s)  _  m&  (s)]®'2* 

. r  n  n  n 

m*0 

where  e'2*  indicates  evaluation  in  the  time  domain  of  the  quantity  in  brackets  at  the  retarded  time 
t-  2. 

Equations  (6.9)  and  (6.12)  constitute  two  equations  for  the  two  unknowns  <l>n(l,s)  and ^(s) 
in  terms  of  the  radial  derivative  <!>.,( l,s).  They  too  are  key  equations  in  the  exact  fluid-structure 
interaction  formulation  of  Section  7. 


(6.11) 


(6.12) 
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6.3  MODAL  DAA  EQUATIONS  FOR  THE  EXTERNAL  FLUID. 


The  definitions  p  =  <j>  and  li  =  -  V<)>  ,  along  with  (6.1)  and  (6.3),  yield  as  the  relation 

between  the  nth  components  of  nondimensional  pressure  and  radial  fluid-particle  displacement 
on  the  unit  sphere 


n  (s)  =  -su  (s)  (6.13) 

kB(Sr» 


where  k„  is  given  in  (6.6).  This  relation  is  the  starting  point  for  the  approximations  derived 
below.  Although  the  derivations  are  carried  out  in  transform  space,  the  resulting  approximations 
are  all  simple  polynomials  in  s,  so  that  inversion  to  the  time  domain,  which  yields  ordinary 
differential  equations,  is  straightforward. 

First,  we  derive  approximations  valid  at  early  time,  which  corresponds  to  s— *»,  and  other 
approximations  valid  at  late  time,  which  corresponds  to  s— >0.  For  s— (6.6)  yields 


k„(s>  -  *  r.,s-  -  r„2s-*  *  ...) 


k.(s)  -  -it’ll  .  r  ,s-'  ♦  (r„,  ♦  r„2)s -  ...] 


(6.14) 


Long  division  then  gives  the  s— »°°  ratio 


k„'(s) 

kjs) 


+  0(s  ’2) 


(6.15) 


For  s — >0,  we  premultiply  (6.6)  through  by  sIMl,  which  yields 


•■"w  -  |=-(r„  *  u  *  r«> 


+  ...) 


(6.16) 


from  which  we  obtain 
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(6.17) 


s-'k.00  -  *  i>r„  *  (r„  *  nr^ps 

Mr.,..,,  *(n  -  i)r.M]s2  *  ...} 


Long  division  then  gives  the  s— >0  ratio 


k„(s) 

K(s) 


-_l_s  +  0(s3) 
n+1 


(6.18) 


First-  and  second-order  early-time  approximations  are  obtained  by  keeping  terms  in  (6.15) 
through  s°  or  s'1,  and  introducing  each  result,  one  at  a  time,  into  (6.13).  The  two  ETA’s  are 

ETAi"(s):  p  (s)  =  s  u  (s) 

(6.19) 

ETA2“(s):  (s  +  l)£n(s)  =  s2un(s) 

Similarly,  first-  and  second-order  late-time  approximations  are  obtained  by  retaining  terms  in 
(6.18)  through  s1  or  s2  and  introducing  each  result,  one  at  a  time,  into  (6.13).  The  results  are 

LTA,“(s):  (n  +  l)p  (S)  -  s*u  (s) 

(6.20) 

LTA2"(s):  same  as  LTA,"(s) 

DAA’s  will  now  be  derived  by  scalar  matching.  As  in  the  previous  sections,  the  method 
consists  of  selecting  a  suitable  trial  equation  in  transform  space  that  contains  an  appropriate 
number  of  arbitrary  coefficients,  and  then  determining  those  coefficients  by  forcing  the  trial 
equation  to  fit  both  the  corresponding  ETA  and  LTA  for  large  and  small  s,  respectively. 

As  will  soon  be  evident,  the  trial  equation  for  the  first-order  DAA  is 

(P,ns  +  P0n)p  (s)  =  s2u  (s)  (6-21) 

■*“n 

where  P,"  and  P0“  are  coefficients  to  be  determined.  This  equation  will  match  the  first  of  (6.19) 
for  arbitrarily  large  s  only  if  P,“  =  1,  and  it  will  match  (6.20)  for  arbitrarily  small  s  only  if 
P0"  =  n  +  1.  Thus  we  have  the  first-order  modal  DAA  in  transform  space 
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DAA,“(s): 


(s  +  n  +  lVp  (s)  =  s2u  (s) 

,4“n 


(6.22) 


and  in  the  time  domain 

DAA,"(t):  p  (t)  +  (n  +  l)p  (t)  =  u  (t)  (6.23) 

A  logical  extension  of  (6.21)  that  satisfies  ETA,"  as  s— and  LTA,"  as  s-40  produces 
the  trial  equation  for  the  second-order  DAA,  viz., 

[s2  +  P,“s  +  (n  +  l)P0n]p  (s)  -  (s3  +  P0Ds2)u  (s)  (6.24) 

**n  ™“b 

To  match  ETA2"(s),  we  divide  this  equation  through  by  s2  and  the  second  of  (6.19)  through  by 
s;  then  we  match  the  resulting  equations  through  order  s'1  to  obtain  P0"  =  P,“  -  1.  To  match 
LTA2"(s),  we  divide  (6.24)  through  by  (n+l)P0“,  divide  (6.20)  through  by  n  +  1,  and  match  the 
resulting  equations  through  order  s;  this  yields  P,"  =  (n+1).  Thus,  introducing  these  results  into 
(6.24),  we  obtain  the  second-order  DAA  in  transform  space  as 

DAA2"(s):  [s2  +  (n  +  l)s  +  n(n  +  1)]£  (s)  ■  (s  +  n)s2u  (s)  (6.25) 

and  in  the  time  domain  as 

DAA2"(t):  £n(t)  +  (n  +  l)£n(t)  +  n  (n  +  l)£n(t)  *  £n(t)  +  nujt)  (6.26) 

6.4  MODAL  DAA  EQUATIONS  FOR  THE  INTERNAL  FLUID. 

For  the  internal  fluid,  the  equation  analogous  to  (6.13)  is 

-iLl £  (s)  =  su  (s)  (6.27) 

i„(s) 

where  i„(s)  is  given  in  (6.11)  and  u(s)  is  defined  for  the  internal  fluid  as  possitive  inward. 
Following  the  procedure  of  Section  6.3,  we  determine  the  singly  asymptotic  approximations  by 
examining  (6.27)  as  s— and  s— »0. 
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For  s-**>,  (6.11)  yields 


i.(s)  -  is-'e‘[l  -  r„s-‘  •  r  ,s  -1 
'.'(»)  -  Is-'e-n  -  (r„,  ♦  i)s-'  -  <r„  ♦  2r.,)s  -J  -  ...] 

Long  division  then  gives  the  s— ratio 

in'(s)  _2, 

i0(s) 


(6.28) 


(6.29) 


For  s— >0,  (6.11)  gives 


where 


m«0 


(6.30) 


[1  -3-5--(2n  +  2m  +  l))2mm! 
Differentiation  followed  by  long  division  then  gives  the  s— >0  ratio 

i  ’(s) 

-1—  *  ns*1  +  2(yn/yn0)s  +  0(s3) 


(6.31) 


(6.32) 


ETA,”  and  ETA2"  are  obtained  by  keeping  terms  through  s°  or  s'1  in  (6.29)  and  introducing 
each  result,  one  at  a  time,  into  (6.27).  The  first  two  ETA’s  are 


ETA,“(s): 

ETA2“(s): 


£  (s)  =  SU  (s) 

— n  n 

(s  -  1) £  (s)  =  S2U  (s) 


(6.33) 


LTA,"  and  LTA2”  are  obtained  by  retaining  terms  in  (6.32)  through  s  or  s2  for  n=0  and  s'1  or  s° 
for  n>0,  and  introducing  each  result,  one  at  a  time,  into  (6.27).  The  results  are 
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LTAj"(s): 


£  (s)  =  3  u  (s) 

o  o 

no  (s)  =  s2u  (s),  n>0 


(6.34) 


Because  the  LTA’s  for  n  =  0  are  distinct  from  their  counterparts  for  n  >  0,  DAA 
derivations  will  be  performed  in  pairs.  The  trial  equation  for  the  first-order  DAA  for  n  =  0  is 

£  (s)  =  (U,°s  +  U0°)u  (s)  (6.35) 

*“0  “O 

where  U,°  and  U0°  are  coefficients  to  be  determined.  This  equation  will  match  the  first  of  (6.33) 
for  arbitrarily  large  s  only  if  U,°  =  1;  it  will  match  the  first  of  (6.34)  for  arbitrarily  small  s  only 
if  U0°  =  3.  Thus  we  have  the  first-order  DAA  for  n=0  in  transform  space 

DAAi“(s)  for  n=0:  £  (s)  =  (s  +  3)u  (s)  (6.36) 

and  in  the  time  domain 

DAA^t)  for  n=0:  £  (t)  *  u  (t)  3  u  (t)  (6.37) 

The  trial  equation  for  the  first-order  DAA  for  n  >  0  is 

(P,“s  +  P0“)£  (s)  =  s2u  (s)  (6.38) 

This  equation  will  match  the  first  of  (6.33)  for  arbitrary  large  s  only  if  P,°  =  1,  and  it  will  match 
the  second  of  (6.34)  for  arbitrary  small  s  only  if  P0”  =  n.  Thus  we  have  the  first-order  DAA  in 
transform  space 

DAA,"(s)  for  n>0:  (s  +  n)£  (s)  =  s2u  (s)  (6.39) 

and  in  the  time  domain 

DAA,n(t)  for  n>0:  £  (t)  +  n£  (t)  =  ii  (t)  (6.40) 

n  *"n  ““o 
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For  the  second-order  DAA,  the  trial  equation  for  n=0  is  constructed  by  extension  of  (6.35) 
as 


(s  +  P0°)p  (s)  =  (s2  +  U,°s  +  3P0°)  u  (s)  (6.41) 

o  ““o 

To  match  ETA2,  we  divide  both  this  equation  and  the  second  of  (6.33)  through  by  s,  and  then 
require  that  the  resulting  equations  match  through  order  s'1.  That  is,  we  require  (1  +  P0°s  '  '^(s) 
=  s(l  +  U,°s '  '^(s)  to  match  (1  -  s '  ‘^(s)  =  su^s)  through  order  s'1;  long  division  and  matching 
yields  U,°  =  Pc°  +  1.  To  match  LTA2,  we  divide  (6.41)  through  by  Po0  and  match  through  order 
s'1  the  resulting  equation  to  the  first  of  (6.34);  this  yields  PQ°  =  2.  With  P0°  and  U,°  thus 
determined,  (6.41)  becomes 

DAA2”(s)  for  n=0:  (s  +  2)p  (s)  =  (s2  +  3s  +  6)u  (s)  (6.42) 

In  the  time  domain,  this  equation  is 

DAA2“(t)  for  n=0:  p  (t)  +  2  p  (t)  =  ii  (t)  +  3  u  (t)  +  6  u  (t)  (6.43) 

Similarly,  by  referring  to  (6.38),  the  trial  equation  for  the  second-order  DAA  for 
n  >  0  is  constructed  as 


(s2  +  P^s  +  nP0“)  p  (s)  =  (s3  +  P0"s2)u  (s)  (6.44) 

n 

By  following  here  the  same  matching  procedure  as  that  applied  to  (6.38),  we  obtain 

P0"  =  Pi”  +  1  and  P,“  =  n.  Thus,  introducing  these  results  into  (6.44),  we  find  the  second-order 

internal  DAA  for  n  >  0  in  transform  space  to  be 

DAA2“(s)  for  n>0:  [s2+  ns  +  n(n  +  1)]d  (s)  *  (s  +  n  +  l)s2u  (s)  (6.45) 

and  in  the  time  domain  to  be 

DAA2”(t)  for  n>0:  p  (t)  +  np  (t)  +  n(n+l)p  (t)  =  u  (t)  +  (n+l)ii  (t)  (6.46) 

n  n  **“n  — n  “n 
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SECTION  7 


TRANSIENT  EXCITATION  OF  A 
FLUID-FILLED,  SUBMERGED  SPHERICAL  SHELL: 
EXACT  AND  DAA  FORMULATIONS 


A  canonical  problem  in  transient  fluid-structure  interaction  is  the  excitation  of  a  submerged 
spherical  shell  by  a  plane  acoustic  wave.  Exact  shell-response  solutions  for  an  empty  spherical 
shell  were  first  provided  by  Huang,  1969;  the  results  produced  by  these  solutions  contained  small 
errors,  which  were  subsequently  found  and  corrected  by  Huang,  et  al.,  1977,  and  by  Geers,  1978. 
In  1979,  Huang  also  obtained  solutions  for  the  problem  of  plane- wave-excited  concentric  spherical 
shells  with  fluid  present  in  the  annular  region  between  the  shells  and  absent  inside  the  inner  shell. 
The  plane-wave  excitation  problem  for  a  submerged  single  shell  filled  with  fluid  has  apparently 
never  been  solved. 

In  this  section,  exact  shell-response  and  acoustic-pressure  solutions  are  obtained  for  a 
fluid-filled,  submerged  spherical  shell  excited  by  a  plane  acoustic  wave.  The  method  of  separation 
of  variables  is  used  to  construct  generalized  Fourier  series  expressions.  For  some  response 
quantities,  especially  surface  pressures,  convergence  of  the  series  is  not  satisfactory,  so  special 
techniques  are  employed,  with  gratifying  results.  As  in  Section  6,  nondimensional  variables  are 
used,  with  length  normalized  to  a,  time  to  a/ce,  and  pressure  to  Pece2,  where  pe  and  ce  pertain  to 

the  external  fluid. 

DAA  Fourier  series  solutions  are  also  obtained  for  the  purpose  of  evaluating  the  doubly 
asymptotic  approximations  for  internal  acoustic  domains  developed  in  Sections  3  and  5. 
Numerical  results  produced  by  the  exact  and  DAA  Fourier  series  are  presented  in  Section  8. 

7.1  DESCRIPTION  OF  THE  PROBLEM. 

A  diagram  of  the  problem  appears  in  Figure  1.  For  generality,  the  internal  and  external 
acoustic  fluids  are  regarded  as  having  different  mass  densities,  pj  and  pe,  and  different  sound 

speeds,  q  and  ce  The  shell  material  is  elastic  and  isotropic,  with  density  Pq  and  plate  velocity  Cq  = 
[E/(l-v2)]1/2,  where  E  is  Young's  modulus  and  v  is  Poisson’s  ratio.  The  shell's  thickness-to- 
radius  ratio  h/a  is  sufficiently  small  that  thin-shell  theory  suffices.  The  response  equations  are 
formulated  in  terms  of  four  field  variables:  the  meridional  and  radial  shell  displacements  v(6,t)  and 
w(0,t),  and  the  internal  and  external  velocity  potentials  <^(r,0,t)  and  <^(r,0,t). 


5  1 


The  expansions  in  generalized  Fourier  series  of  the  internal  and  external  velocity  potentials 
are  given  by  (6.1),  and  those  of  the  shell  displacements  are  given  by 

oo 

v(0,t)  =  -  I  Vn(t)^Pn(COS0) 
n-1 

(7.1) 

oo 

w(0,t)=  X  Wn(t)Pn(COS0) 
n»0 


Acoustic  pressure  and  radial  fluid-particle  velocity  at  the  external  shell  surface  are  related  to  the 

e  e  ,e  e 

external  velocity  potential  there  by  p  =  <f>  and  ji  *  -  <£  ,  where  a  dot  denotes  partial 

~  ,r 

differentiation  with  respect  to  nondimensional  time  and  the  r-subscript  denotes  partial 
differentiation  with  respect  to  the  nondimensional  radial  coordinate;  we  recall  that  an  underline 

means  evaluation  at  r  =  1.  The  corresponding  relations  for  the  internal  acoustic  medium  are 

i  i  i  i  t  • 

p  =  (p/p^  and  ji  =  4>  ,  where  u1  is  positive  inward. 

,r 

7.2  MODAL  EQUATIONS  OF  MOTION  FOR  THE  SPHERICAL  SHELL. 

The  equations  of  motion  for  the  nth  Fourier  component  of  shell  response  may  be  written 
[Junger  and  Feit,  1972] 


°(n  +  l)vn  +  An  vn  +  An  wn  =  0 

VW  U/W 

*n+  *n  %  +  *n  wn  =  MPn 


(7.2) 


in  which  p  =  (pe/po)(a/h),  pn  is  the  nth  Fourier  component  of  net  pressure  acting  radially  outward 


on  the  shell,  and 

W 

An  =  n(n  +  1)(1+  e)5nY0 

VW 

An  =  n(n  +  1)(1+  v+  e$n)y0 

WW 

An  =  [2(1+ v) +n(n  +  l)e5n]Y0 


(7.3) 


where  e  =  (h/a)2/12,  yQ  =  {c^c^}2,  $n  =  n(n+l)-l+v. 
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7.3  EXACT  FLUID- STRUCTURE-INTERACTION  EQUATIONS. 


Force  compatibility  at  the  surface  of  the  shell  requires  that  p  =  p1  -  pe  ,  where  p1  and 

n  n  n  n 

pe  are  the  pressures  exerted  by  the  internal  and  external  fluid,  respectively.  But  pe  is  the  sum  of 
n  n 

the  known  incident  pressure  p°  and  an  unknown  scattered  pressure  ps  ;  furthermore,  modal 

n  n 

surface  pressures  are  related  to  corresponding  velocity  potentials  as  p1  =  (p/pjtft^ 

n  1 

e  «e 

and  p  =  <!>_.  Hence  p  may  be  expressed  as 
n  ~ n 


P  “(P/Pe)^1  ~  p°  ■  $ 


s 

> 

'  n 


(7.4) 


With  modal  radial  fluid-particle  velocities  related  to  corresponding  velocity  potentials  as 

. i  i  J  .S  ,S  .e  .o  .s 

a  =  <p  and  a  *  -  <p  ,  and  with  u  =  u  +  u  ,  geometric  compatibility  at  the  inner  and 

n  u  ^  Q  OyT  non 

outer  surfaces  of  the  shell  requires  that 


w 


n 


nj 


(7.5) 


Note  that  circumferential  geometric  compatibility  is  not  enforced,  as  the  fluid  is  inviscid. 

Now  (6.9)  for  the  internal  fluid  and  (6.4)  for  the  external  fluid  enable  us  to  eliminate  $ 

n,r 

s 

and  $  from  (7.5),  which  yields 

n.r 

(cJcW  -4?1  -  V1  =  -  wn 

1  n  n  n 

(7.6) 

,s  S  s  . 

+§  +  lj5  =Wa-U 
nan  n 


s 

Also,  we  observe  from  (6.7)  that  each  scattered  residual  potential  !£  (u  is  related  to  the 

n 

s 

corresponding  velocity  potential  $  (t)  through  the  differential  equation 


I  rnmSS 


m-0 


n,n-m 


m-1 


s 

n  ,n  -m 


(7.7) 
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where  the  subscript  n-m  denotes  (n-m)-fold  differentiation  in  time.  Finally,  we  see  from  (6.12) 
that  each  internal  residual  potential  4^(0  is  related  to  the  corresponding  velocity  potential  (0 
and  historical  data  through  the  delayed-differential  equation 


1  <-  -  I  <  -  Dmm(cyc  i)"'"rim4‘ 

*  n  n-  m  *  r 


n,n- m 


n  ,n  -m 


n  _ 

+  (  -  1)“  L  (c(/c.)n”mrB111  %  _  -  2(c€/ci)<j>1  -m^1 

m„0  L  n,n-m  n.n-m+l 


■'  1 
n  ,n-mj 


t-2c./c 

•  t 


7.4  ASSEMBLY  OF  THE  EXACT  RESPONSE  EQUATIONS. 

The  ensemble  of  equations  (7.2),  (7.4),  (7.6),  (7.7)  and  (7.8)  constitutes  a  set  of  seven 

s  s  i  i 

equations  for  the  seven  unknowns  v_(t),  w_(t),  p  (t),  <J>  (t),  ig  (t),  <j>  (t)  and  (0, 

n  n  n  n  n 

given  the  incident-wave  functions  p°(t)  and  li  (t) .  However,  (7.8)  presents  a  problem  in  that 

„  n 

the  highest  derivative  of  ^  (t)  in  the  time-delayed  term  would  be  one  order  higher  than  the 

n 

highest  derivative  of  4>  (t)  previously  calculated.  This  problem  can  be  overcome  by  numerical 

n 

differentiation  of  the  nth  derivative  of  4>  (t),  which  is  not  particularly  appealing.  Fortunately,  it 

n 

can  also  be  overcome  by  (n-m)-fold  differentiation  of  the  first  of  (7.6),  which  permits  the 
i 

replacement  of  4>  ^  in  (7.8)  by  means  of  the  relation 


i  i  i 

(c  Jc  .)<b  =  4>  +  V  ~  w 

n,n-m  n  ,n-m  n,n-m+i 


Unfortunately,  this  equation  brings  with  it,  for  n  >  0,  unacceptably  high  derivatives  of  wn,  which 

S  s  i 

we  avoid  by  introducing  the  integrated  variables  Vn(t),  Wn(t),  $  (t),  3?  (t) ,  $  (t)  and 


3?  (t),  defined  by 
n 
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(7.10) 


v  =  V 
v  n  T  n 


,n  wn  ^n,n  ~ 


c  s  11  ;  1 

mS  =  S'  <t)  =  Cp  m  =% 

*n  n.n  *n  n,n 


n,n 


Thus,  inserting  (7.4)  into  (7.2),  integrating  (7.7)  and  (7.8)  n  times,  and  introducing  the  integrated 
variables  into  the  resulting  equations  and  into  (7.6),  we  obtain  the  following  six  equations  for  six 
unknowns: 


VW 


Vn.n+  *»■“ 


W„.„*2  +  Om  +  »X  *  MliSn  n+1  -  (P/P^‘n  ,,Hl  -  -  MP 


w  ..-4s  -s‘  -i° 

n,n+l  n,n+l  n.n  n,n  n 


W 


+  (c7c.)$‘  -V  =0 

n.n+l  e  i  nai+i  n,n  n,n 


(7.11) 


£  rnm(m^s  -  n?s  )  =  0 

*-*  nm'  n.n-  m  n,n-m 
m=0 


y  (-  ^(c./c/rjm^1  -X  )* 

'  i  &  nm'  n,n-m  n,n- m 

n  r  .  , 

(  -  1)”  I  (c./cJ^1^nmhi!,  +(m  +  2)^‘  -2W 

'  7  **  '  i  &  nm[^  n,n- m  n,n- m  n,n-m+lj 


m“0 


m-0 


t-2c./c, 


Once  the  solutions  to  these  equations  have  been  obtained  for  several  values  of  n  and  the  desired 
modal  response  histories  have  been  determined  in  accordance  with  (7.10),  Fourier  superposition 
yields  response  histories  for  surface  pressures  and  shell  responses  at  desired  locations  in 
accordance  with  (6.1)  and  (7.1). 
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7.5  ASSEMBLY  OF  THE  DAAj  RESPONSE  EQUATIONS. 


To  obtain  a  set  of  equations  corresponding  to  (7.1 1)  that  are  based  on  a  DAAj  treatment  of 
the  fluid-structure  interaction,  we  first  use  (6.37)  and  (6.40)  for  the  internal  fluid  along  with  the 
compatibility  relation  =  -wn  to  get  [cf.  the  first  of  (7.6)] 


(c«/c  ^  -  -  [w  o  +  3(0  70  e)' w0] 

(c^/cV  +  n<j/  =  -  wn  (n  >  0) 
1  n  n 


(7.12) 


•  O  s  • 

Next,  we  use  (6.23)  for  the  external  fluid  along  with  the  compatibility  relation  un  +  un  =  wn  to 
obtain  [cf  the  second  of  (7.6)] 


s  s  .  .o 

<p  +  (n  +  l)d>  =  w  -  _u 
~  n  n  n 

Finally,  we  introduce  (7.4)  into  (7.2)  to  eliminate  p  as  an  unknown. 

n 

s  i 

equations  for  the  four  unknowns  vn(t),  w_(t),  ?  (0  and  4>  (t): 

n  n 

w  vw 

n(n  +  l)vn  +  Xnvn+  Xa  w„  =  0 

VW  WW  S  1 

^n+An  vn  +  *n  wn  +  ~  (P/P^  J=“PP 

s  s  .o 

w  -  $  -  (n  +  1)$  =  u 
n  n  n 


(7.13) 

Thus  we  obtain  four 


O 

n 


(7.14) 


w0+  (cfc  )$'  +  3(c./ce)wo  =  0 

1  O  1 

wn+(ce/c  )^  +  nc]>  =0  (n  >  0) 

1  n  n 
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7.6  ASSEMBLY  OF  THE  DAA2  RESPONSE  EQUATIONS. 


To  obtain  a  set  of  equations  corresponding  to  (7. 1 1)  that  are  based  on  a  DAA2  treatment  of 
the  fluid-structure  interaction,  we  first  use  (6.43)  and  (6.46)  for  the  interior  fluid  along  with  the 
compatibility  relation  =  -wn  to  get  [cf.  (7.12)  and  the  first  of  (7.6)] 

(cJc.W  +  24/  =  -  [w0 +  3(c./ce)w0  +  6(c./ce)2w0] 

1  o  o  1  1 

(cJcW  +n4>‘  +n(n  +  -  -  [wn  +  (n  +  lHc./c^wJ 

*  n  n  *11  1 

Next,  we  use  (6.27)  for  the  external  fluid  along  with  the  compatibility  relation  =  wn  to 

obtain  [cf.  (7.13)  and  the  second  of  (7.6)] 

<J>8  +  (n  +  1)4>S  +  n(n  +  1)<|>S  -  w_  -  i°  +  n(w  -  i°) 

-n  -n  -n  “  n  ”  n  (7.16) 

Finally,  we  introduce  (7.4)  into  (7.2)  to  eliminate  p  as  an  unknown.  This  yields  the  following 

n 

s  i 

four  equations  for  the  four  unknowns  vn(t),  wn(t),  SM*)  and  <J>^(t): 

vv  vw 

n(n  +  l)vn  +  *n  v„  +  *n  wn  =  0 

VW  ww  s  1  0 

wn  +  Xn  Vn  +  Xn  wn  +  ^  ”  (P/Pjg  1  m  "  PP 

n  *  n  n 

s  s  s  „o  .0 

w  -  ®  +  nwn  -  (n  +  1)4>  -  n(n  +  1)4>  =  +  nu 

n  “n  n  “n  "n  “  ®  (7.17) 

w0+(cAi)^  +  3(c  /ce)w0+2^'  +  6(c /ce)2wo-0 

1  O  1  O  1 

wn  +  (c  Jc . )$'  +  (n  +  lXc./Cg)^  +  +  n(n  +  lKc./c^  *0  (n  >  0) 

In  contrast  to  (7.1 1),  (7.14)  and  (7.17)  are  low-order  ordinary  differential  equations  in  the 
direct,  not  integrated,  variables.  Once  modal  solutions  to  these  equations  have  been  obtained  for 


(7.15) 

(n>l) 
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several  values  of  n,  Fourier  superposition  yields  results  for  surface  pressures  and  shell  responses 
at  desired  locations  in  accordance  with  (6.1)  and  (7.1). 

7.7  MODIFIED  CESARO  SUMMATION  FOR  IMPROVED  CONVERGENCE. 


As  mentioned  above,  the  generalized  Fourier  series  calculated  for  some  of  the  response 
quantities  of  interest  do  not  converge  satisfactorily.  This  is  certainly  to  be  expected  in  response 
histories  that  contain  discontinuities,  where  pronounced  non-physical  oscillations  appear  (Gibb’s 
phenomenon).  A  superposition  technique  that  has  proven  effective  in  reducing  these  oscillations  is 
due  to  Cesaro  (Apostol,  1957). 

With  Sn  as  the  partial  sum  of  an  infinite  series  through  the  first  N+l  terms,  the  Nth  Cesaro 
sum,  0^,  of  that  series  is  the  arithmetic  mean  of  the  first  N+l  partial  sums  Sj^,  i.e.. 


N 


1 

N  +  1 


N 


ISM 

M-0 


(7.18) 


Introducing  the  first  of  these  into  the  latter  and  expanding,  we  find  that  the  Cesaro  sum  may  be 
written  explicitly  as 


fT  —  v  +  -  —  Y  *  Y  +  .  •  .  4.  A  Y 

N  x0  N  +  1  X1  N+l  2  N  +  l  N  (7.19) 

A  useful  interpretation  of  partial  and  Cesaro  summation  consists  of  regarding  each  as  a 
'  gital  weighting  filter  for  an  infinite  series.  In  this  interpretation,  partial  summation  employs  unit 
weights  for  the  first  n+l  terms  and  zero  weights  for  the  rest,  and  Cesaro  summation  employs 
weights  that  decrease  linearly  from  one  for  the  first  term  to  zero  for  xn+i  and  beyond.  Clearly,  the 

filter  characteristic  for  Cesaro  summation  is  more  gradual  than  that  for  partial  summation. 

Standard  Cesaro  summation  is  not  appropriate  in  the  present  problem  because  weights  less 
than  unity  for  n=l  and  n=2  produce  inaccurate  late-time  asymptotic  results  for  translational  velocity 
and  deformational  displacement  of  the  shell  [Geers,  1974].  The  procedure  is  easily  modified, 
however,  to  produce  unit  weighting  for  modes  0,  1 ,  and  2,  and  linearly  decreasing  weights  for 
modes  3  through  N.  The  resulting  filter  characteristics  for  N  =  5  and  N  =  8  are  shown  in  Figure 
2;  also  shown  are  the  corresponding  filter  characteristics  for  partial  summation. 

Cesaro  summation  can  dramatically  improve  convergence,  as  demonstrated  in  Figures  3 
and  4.  The  figures  show  pressure  histories  generated  for  a  free-field  step- wave  by  the 
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superposition  of  modal  pressure  histories  in  accordance  with  (6.1).  Each  modal  pressure  history 
is  given,  for  0  <  t  <  2,  by 


JT 

2  (t)  =  (n  +  ^)PjjH(t  -  cos  6  -  1)  Pn(cos0)  sin  0d0 

0  (7.20) 

For  t  >  2,  p  ( t)  =  p  (2) .  The  integral  in  (7.20)  is  easily  evaluated  by  the  change  of  variable  t  = 
~n  — n 

cos0.  Pressure  histories  are  shown  at  three  points  on  the  locus  r  =  1  in  the  spherical  geometry. 
The  true  histories,  of  course,  are  step-functions,  with  discontinuities  at  t  =  0,  1  and  2  for  6  -  7r, 
rt/2  and  0,  respectively. 

Figures  3  and  4  also  provide  values  of  integrated  mean-square  error,  given  by 


2 


where  P  (t)  is  the  summed  history  and  P  (0  is  the  exact  history.  The  values  indicate  that. 
In  Ex 

while  modified  Cesaro  sums  may  be  superior  to  standard  partial  sums  at  some  points,  they  may  be 
inferior  at  others.  This  is  demonstrated  more  comprehensively  in  Figure  5,  which  shows 
integrated  mean-square  error  characterizing  modified  Cesaro  and  standard  partial  sums  for  step- 
wave  pressure  histories  on  the  locus  r  =  1.  As  one  would  expect  from  the  overall  optimality 
attribute  of  Fourier  series,  the  average  integrated  mean-square  error  for  standard  partial  summation 
is  less  than  that  for  modified  Cesaro  summation.  However,  the  maximum  error  produced  by  the 
latter  is  less  than  that  produced  by  the  former.  Furthermore,  standard  partial  summation  incurs  its 
largest  errors  at  e  =  0  and  8  =  n,  which  are  often  the  points  of  greatest  interest. 

From  the  preceding,  modified  Cesaro  summation  yields  for  shell  velocities  and  surface 
pressures 

N 

v(0,t)  =  '  I  Cnvn(t)-jjPn(cos0) 

n  *1 

N 

w(0,t)  =  £  CnWn(t)Pn(COS0) 

n-0 

(7.22) 
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where  CQ  =  C j  =  C2  =  1  and  Cn  =  (N+l-n)/(N-l)  for  n  >  2. 


7.8  PARTIAL  CLOSED-FORM  SOLUTION  FOR  IMPROVED  CONVERGENCE 


It  is  clear  from  Figures  3  and  4  that  no  superposition  of  modal  solutions  can  reproduce  the 
jump  in  pressure  at  a  discontinuous  wave  front  This  deficiency  is  even  more  pronounced  in  the 
vicinity  of  the  point  of  first  contact  between  the  wave  front  and  the  spherical  shell,  where  the 
pressure  initially  doubles.  Hence  we  introduce  here  a  method  to  alleviate  this  convergence  problem 
by  assembling  the  complete  solution  as  the  sum  of  a  closed-form  initial  solution  and  a 
complementary  series  solution. 

The  first  step  in  the  method  involves  retention  of  the  terms  in  (7.1 1)  that  dominate  early- 
time  response;  this  yields  the  initial-response  equations 
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(7.23) 


W*  a..  +  (c7c.)<t>1  =0 

n.n+l  e  njl+1 


where  the  asterisk  denotes  initial-response  quantities.  It  is  apparent  that  these  equations  neglect  all 
stiffness  effects  in  the  shell  and  invoke  the  first-order  early-time  (plane-wave)  approximation  for 
the  fluid-structure  interaction.  Because  they  do  not  involve  the  modal  index  n  as  an  explicit 
parameter,  they  can  be  summed  in  closed  form.  Thus,  eliminating  the  two  integrated  velocity 
potentials  from  the  first  equation  by  using  the  other  two,  and  then  utilizing  (7.10),  we  may  write 

s*  i* 

the  summed  equations  in  toms  of  the  direct  variables  w*(0,t),  $  ( 0,  t)  and  $  ( 0,  t)  as 
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w*  +  uw*=  -  (j(2°  -  u^J 


5  =  w*  -  _u 


(7.24) 


$  “-(c./c^w* 

where  o  =  |u(  1  +  piq/peeg). 

The  second  step  consists  of  solving  (7.24)  for  a  prescribed  incident  wave.  For  example, 
an  incident  step  wave  propagating  to  the  right  that  at  t  =  0  first  contacts  the  shell  at  6  -  7r  yields 


j>°(0,t)  *  PjH(t  -  cos  0-1) 
jj  (0,t)  =  PjCOS0H(t  -  cos0  -  1) 


(7.25) 


where  H( )  is  the  Heaviside  step-function.  For  this  excitation,  the  closed-form  solutions  to  (7.24) 
with  quiescent  initial  conditions  are 

^*(0,t)  =  -mu-1Pt(1  -cos0)[l  -e~u(t~cos0  I)]H(t  -cos0  -  1) 


£S*(0,t)  =  -  piT^l  -e'u(t_C08e“1)]H(t  -cos0  -  l) 

-  P  cos0{l  -  pu  _1[1  -  e‘u(1_cos0"!)]}H(t  -  cos0  -  1) 


(7.26) 


<j)**(0,  t)  *  pu _1(c  ./c^Pjfl  -  cos  0)[1  -  e"v(t'cos0'1)]H(t  -  cos 0  -  1) 

The  third  step  requires  that  the  closed-form  solutions  be  used  to  compute  modal  initial- 
response  histories  from  the  standard  formula 


=  (n  +  ^)J{w*(0,t),<j>S(0,t),  <j>  (0,t)  }Pn(cos0)  sin0  d0 


i  e .  *  s*  »' 

*  \  I  f  •  _  /A  A.  \  /  A  .  \  1 


(7.27) 
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Then  these  histories  are  numerically  integrated  in  time  to  tabulate  the  modal  initial  responses 

w!(t),  £S*(t ),$S*  (t),-  •  •,  $S*(t),  4‘(t),  ^‘(tXi1*  ox-  •  \  £s*(tx 

n  n,n  n,n-l  n,l  n  n’n  n,n-l  n,l  n 

* 

The  fourth  step  involves  associating  with  each  modal  initial  response  [e.g.,  wn(t)J  a  modal 

complementary  response  \e.g.,  w*(t)]  such  that  the  sum  of  the  two  yields  the  true  modal  response 
*  + 

[e.g.,  wn(t)  =  wn(t)  +  wn(t)],  subtracting  each  of  (7.23)  from  its  counterpart  in  (7.1 1),  expressing 

the  last  two  of  (7. 1 1)  in  terms  of  initial  and  complementary  responses,  and  invoking  the  second  of 
(7. 10).  This  yields,  for  t  <  2cg/q.  the  complementary-response  equations 
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(7.28) 
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The  final  step  consists  of  the  addition  of  initial  and  complementary  solutions,  and  the  use  of 
modified  Cesaro  summation,  which  yields  for  shell  velocities  and  surface  pressures  [cf.  (7.22)] 

oo 

V(6,t)  -  -  £  Cnvn(t)-|Pn(cose) 


w(0,t)  =  w*(0,t)  +  ICnwn(t)Pn(cos0) 


(7.29) 
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where,  again,  CG  =  C}  =  C2  =  1  and  Cn  =  (N+l-n)/(N-l)  for  n  >  2.  Because  the  delayed  term  in 
the  last  of  (7.11)  was  dropped  in  the  process  of  deriving  (7.28),  (7.29)  are  only  valid  for  t  < 

2cg/cj.  This  is  generally  satisfactory,  in  that  ps(0,t)  contains  no  discontinuity  for  0  >  n/2 

[Friedlander,  1958],  which  the  incident  wave  front  reaches  at  t  =  1 .  Hence  it  is  appropriate  to  use 
(7.29)  for  t  less  than  unity  or  2cg/q,  whichever  is  smaller. 

The  first-order  early-time  (plane-wave)  approximation,  on  which  (7.29)  are  based,  is  only 
accurate  for  t «  1.  The  region  of  validity  of  this  approximation  is  readily  assessed  by  noting  that 
the  second  of  (7.26)  predicts  a  wave-front  jump  in  scattered  pressure  of  -Pjcos0o  at  the  circle  on 

the  shell  defined  by  r  =  1, 0  =  arccos  (t  - 1).  In  contrast,  at  points  on  the  shell  reached  first  by  the 

incident  wave,  the  true  jump  is  unity  [Geers,  1972].  Hence  we  would  expect  (7.29)  to  predict 

discontinuous  scattered- wave  response  accurately  over  the  region  180°<  0  <  155°. 

Partial-closed-form  solution  with  modified  Cesaro  summation  is  also  used  to  obtain  DAA 

* 

results.  The  initial-response  solutions  for  DAAj  are  (7.26),  but  only  wR(t)  needs  to  be  tabulated 
using  (7.27).  The  complementary-response  equations  are 
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The  initial-response  solutions  for  DAA2  are  (7.26)  also,  but  the  only  modal  initial  responses  that 

* 

must  be  tabulated  using  (7.27)  are  wfl(t)  and  wn(t).  The  complementary-response  equations  are 
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+  (c JCy)$l  Mn+lKciCgJw/n^1  +  n(n  +  lKc./c^' 


=  -(c./ce)[wn-n(n  +  lMc./c^w*]  (n  >  0) 


The  partial-closed-form  DAA  solutions  are  obtained  in  accordance  with  (7.29)  for  t  <  1. 
7.9  INTERNAL  ACOUSTIC  FIELDS. 


From  (6.8)  and  from  p*(r,0,t)  =  (pj/pe)4^(r,0,t)  and  u*(r,0,t)  =  -V^V^t),  the  Fourier 

components  of  Laplace-transformed  acoustic  pressure,  radial  fluid-particle  velocity,  and  meridional 
fluid-particle  velocity  at  any  point  in  the  internal  fluid  may  be  expressed  as 


p|,(r,s)  =(p./pe)sf1n(s)iD(rsc^c.) 
su‘n(r,s)  =  -  (sc^c.)f,n(s)in(rsc^c.) 


(7.32) 


sulen(Ls)  =  -  r_1  fin(s)in(rsce/c.) 


where,  again,  in( )  is  the  nth-order  modified  spherical  Bessel  function  of  the  first  kind  and  i^  ( ) 
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is  the  derivative  of  in( )  with  respect  to  its  argument.  But  the  rah  Fourier  component  of  inward 
normal  surface  displacement  is  just  ^(s)  =  -  uj^(l,s).  Hence  we  may  use  (7.32)  to  relate  each 

Fourier  component  of  the  internal  acoustic  field  to  the  corresponding  component  of  inward  normal 
surface  displacement  as 


i  ( rsc  J  c  ) 

P°(r’S)  -  (Pic  i/Pec  e>  <(s) 
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(7.33) 


The  inverse  transform  of  each  of  these  relations  is  a  delayed-differential  equation,  which  can  be 
solved  numerically.  With  solutions  for  several  values  of  n  thus  obtained,  Fourier  superposition 
yields  response  histories  for  the  internal  fields  as 

00 

pl(r,e,t)=  X  Pn(rit)Pn(cos8) 
n-0 


u‘(r,0,t)=  £  u‘n(r,t)Pn(cos0) 

n-0  (7.34) 


Ug(r,0,t)  =  X  uin(r’t)^-Pn(cos0) 

n«l 


It  is  particularly  easy  to  obtain  response  histories  at  the  center  of  the  internal  fluid 
domain, viz.,  at  r  =  0.  At  this  point,  symmetry  arguments  allow  us  to  conclude  that  only  the  n  =  0 
component  of  surface  displacement  contributes  to  pressure  and  only  the  n  =  1  component 
contributes  to  fluid-particle  displacement.  For  these  modes,  inverse  transformation  of  the  first  and 
second  of  (7.33)  yields,  with  geometric  compatibility  requiring  that  ^(t)  =  -  wn(t)  and  with 

integrated  variables  avoiding  high  derivatives, 
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(7.35) 


4>'o(0,t)  -(c./ce)<J)lo(0,t)=  -2w0(t  -C^C.)  ~ 4>o(0,t  -  cjc.) 


Url(0,t)  -  2(c./ce)U1rl(0,t)  +  2(ci/ce)2u‘1(0,t)  =  f-(c</c.)w1(t  -cjc.) 


+  Ufl(0,t  -2Ce/c.)  +  2(c./Ce)Ufl(0,t  -2cjc.)  +2(c./Ce)  Url(0,t  -  2c Jc  .) 

where  U*  j  (0,t)  is  defined  by  (jJ.j(0,t)  =  uj.j  (0,t).  From  the  solutions  to  these  equations,  we  obtain 
for  pressure  and  fluid-particle  velocity  at  r  =  0 

p‘(0,0,t)  «  (p./p^fO.t) 


ur(0,0,t)  =  Url(at) 


(7.36) 


u0(O,O,t)  =0 
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SECTION  8 


NUMERICAL  RESULTS 

Computed  transient  response  histories  are  presented  in  this  section  for  a  steel  spherical 
shell  filled  with  and  submerged  in  water,  and  excited  by  a  plane  step- wave.  The  thickness-to 
radius  ratio  is  h/a  =  0.01,  the  mass-density  ratio  is  po/p  =  7.7,  and  the  sound-velocity  ratio  is  cjc 
=  (13.8)1/2.  Fifth-order  Runge-Kutta  integration  was  used  to  solve  (7.11),  (7.14),  (7.17),  (7.28), 
(7.30),  (7.31)  and  (7.35),  and  Simpson’s  rule  was  used  to  perform  the  numerical  integrations  in 
space  and  time  discussed  in  connection  with  (7.27).  Partial  closed-form  solution  with  modified 
Ces^ro  summation,  i.e.,  (7.29),  was  used  during  0  <  t  <  2,  and  modal  superposition  with  modified 
Ces^ro  summation,  i.e.,  (7.22),  was  used  during  2  £  t  <  10,  As  many  as  nine  Fourier  meriodonal 
harmonics  were  used;  response  histories  for  both  N  =  5  and  N  =  8  are  displayed  in  the  figures 
to  indicate  the  degree  of  modal  convergence. 

The  computer  coding  was  validated  by  first  performing  response  computations  for  a  shell 
with  mass-density  and  sound-velocity  ratios  of  unity.  The  results  showed  the  shell  to  be 
essentially  transparent,  as  desired.  Computations  were  also  performed  for  an  empty  submerged 
shell,  and  the  results  were  compared  with  the  previous  results  of  Geers,  1978;  excellent 
agreement  was  found.  Finally,  early-  and  late-time  responses  were  checked  against  easily 
calculated  early-time  jump  and  late-time  static  values;  proper  behavior  was  obtained  in  all  cases. 

8.1  EXACT  RESULTS. 

Response  histories  corresponding  to  an  exact  treatment  of  the  fluid-structure  interaction 
are  displayed  in  Figure  6  through  15.  Figure  6  shows  surface-pressure  histories  at  0  =  n  that 
were  generated  without  the  use  of  partial  closed-form  solution;  the  absence  of  the  requisite  jump 
in  and  the  poor  convergence  for  0  £  t  <  2  are  clearly  seen.  Figure  7  shows  the  same  histories 
generated  with  partial  closure;  the  improvement  is  obvious. 

The  pressure  histories  of  Figure  7  exhibit  considerable  texture  produced  by  complex  wave- 
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propagation  effects.  The  external-pressure  history  drops  rapidly  following  its  initial  jump  to  two, 
which  is  characteristic  of  a  thin  shell  (Geers,  1972);  at  late  time,  approaches  its  late-time 
asymptote  of  unity.  The  internal-pressure  history  initially  rises  rapidly  from  zero,  almost 
reaching  the  external-pressure  history,  and  then  oscillates  irregularly  about  its  late-time  asymptote 
of  0.52. 

Figure  8  shows  surface-pressure  histories  at  the  deep-shadow  point  0  =  0.  Virtually  no 
response  is  seen  until  t  =  0.75,  which  demonstrates  the  effectiveness  of  Ces&ro  summation  and 
partial  closure  for  modal  convergence  enhancement.  After  the  arrival  of  stress  waves  in  the  shell, 
which  generate  shell  motion  and  thus  surface  pressures  and  before  the  arrival  at  t  =  2  of  the 
acoustic  wave  that  travels  through  the  internal  fluid,  the  external-pressure  history  exhibits  a 
modest  hump  and  the  internal-pressure  history  shows  a  substantial  negative  excursion.  At  t  = 
2,  both  pressure  histories  rise  rapidly  and  then  oscillate  irregularly  about  their  late-time 
asymptotes  of  1.0  and  0.52,  respectively. 

Shell  radial-velocity  histories  at  8  =  it  and  0  =  0  are  displayed  in  Figure  9.  In  keeping 
with  the  plane- wave  approximation,  the  latter  history  matches  the  internal  surface-pressure  history 
in  Figure  7  at  early  time;  the  two  are  even  quite  similar  until  t  =  2.  A  brief,  but  prominent, 
oscillation  appears  in  the  0  =  7t  velocity  history  during  2.5  <  t  <  4,  which  is  caused  by  the 
internal  pressure  wave  reflected  back  from  the  rear  portion  of  the  shell.  The  motion  after  t  =  4 
consists  of  low-frequency  oscillation  about  the  asymptotic  value  of  0.87  (the  filled  shell  is 
slightly  negatively  buoyant).  As  expected,  no  response  is  seen  at  0  =  0  until  t  =  0.75,  when  shell 
stress  waves  arrive.  The  velocity  history  shows  a  modest  hump  during  0.75  £  t  <  2  (c/.  the  & 
history  in  Figure  8),  after  which  it  rises  rapidly  to  a  peak  value  at  t  =>  3  that  exceeds  the  incident 
wave’s  fluid-particle  velocity  by  nearly  50%.  Subsequent  motion  consists  of  low-frequency 
oscillation  about  the  0.87  asymptote. 

It  is  interesting  to  compare  transient  response  histories  for  a  fluid-filled  shell  with  their 
counterparts  for  an  empty  shell,  as  done  in  Figure  10  through  14.  External-surface  pressure 
histories  at  0  =  7t,  which  are  displayed  in  Figure  10,  exhibit  rather  modest  differences,  with  the 
empty-shell  history  initially  dropping  from  two  at  twice  the  rate  of  the  filled-shell  history, 
reaching  a  substantially  lower  minimum,  and  exhibiting  more  oscillation  at  late  time.  Differences 
are  also  modest  at  0  =  rc/2.  as  seen  in  Figure  11;  here,  the  filled-shell  history  rises  rapidly  after 
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t  =  1,  which  the  empty-shell  history  does  not  do,  and  the  latter  exhibits  more  oscillation  at  late 
time.  Substantial  differences  between  empty-shell  and  filled-shell  pressure  histories  appear  at  0 
=  0,  however,  as  shown  in  Figure  12.  The  empty-shell  history  rises  rapidly  at  t  *  0.75;  the  filled- 
shell  history  rises  even  more  rapidly  at  t  =  2.  Also,  the  empty-shell  history  exhibits  large 
oscillation  at  late  time,  which  is  not  matched  by  filled-shell  history.  In  all  three  figures,  the 
empty-shell  histories  are  quite  smooth  after  t  =  3.5,  while  the  filled-shell  histories  exhibit 
continuing  wave  reflections  in  the  internal  fluid;  in  all  cases,  of  course,  the  late-time  asymptote 
is  unity. 

Velocity  histories  for  the  fluid-filled  and  empty  shells  are  compared  in  Figure  13  and  14. 
In  the  former  figure,  which  pertains  to  0  =  7t,  the  empty-shell  history  initially  rises  at  twice  the 
rate  of  its  filled-shell  counterpart,  reaches  a  substantially  higher  initial  peak,  and  then  oscillates 
with  much  higher  amplitude  about  an  asymptotic  value  of  2.05  (the  empty  shell  is  very  positively 
buoyant).  At  0  =  0,  to  which  Figure  14  pertains,  the  empty-shell  history  rises  immediately  after 
t  *  0.75  to  reach  peak  values  more  than  triple  the  fluid-particle  velocity  of  the  incident  wave. 
Both  figures  fulfill  the  expectation  that,  at  late  time,  the  empty  shell  oscillates  at  a  frequency 
higher  than  that  characterizing  the  filled-shell  oscillation;  also,  the  empty-shell  histories  are 
smoother  than  the  filled-shell  histories  after  t  =  2. 

Finally,  Figure  15  shows  acoustic-pressure  and  fluid-particle- velocity  histories  at  the 
center  point  in  the  internal  fluid.  Plane- wave  propagation  governs  for  1  <  t  <  1.3,  during  which 
time  the  pressure  and  velocity  histories  coincide,  but  the  two  histories  diverge  rapidly  after  that. 
At  no  time  does  the  shell  appear  transparent  to  the  incident  wave,  the  pressure  and  velocity 
histories  for  which  would  appear  in  the  figure  as  a  unit  step-function  at  t  =  1.  Strong  focusing 
effects  at  t  =  3,  5,  7  and  9  are  apparent  in  the  figure,  and  the  two  histories  properly  approach 
their  late-time  asymptotes  of  0.87  and  0.52. 

8.2  DAA  RESULTS. 

Response  histories  corresponding  to  DAA  treatments  of  the  fluid-structure  interaction  are 
exhibited  in  Figure  16  through  24  .  At  0  =  n,  to  which  Figure  16  pertains,  DAA  results  for 
external-surface  pressure  are  quite  satisfactory,  with  the  DAA2  history  following  the  exact  history 


69 


very  closely  everywhere  except  for  2.5  <  t  <  4,  during  which  the  exact  history  exhibits  its  brief 
oscillation.  As  seen  in  Figure  17,  DAA  performance  is  not  as  good  for  internal  pressure,  with 
2  £  t  £  4  being  the  period  of  greatest  difficulty,  also  because  of  the  brief  oscillation.  Here  DAA2 
clearly  emerges  as  superior  to  DAA„  reaching  the  initial  peak  with  greater  accuracy  and 
exhibiting  the  correct  frequency  of  oscillation  at  late  time;  DAA2  does  not,  however,  capture  the 
fine  texture  of  the  exact  history. 

External-  and  internal-surface  pressure  histories  at  0  =  n/2  are  shown  in  Figures  1 8  and 
19.  DAA  performance  is  quite  satisfactory,  with  DAA2  doing  slightly  better  than  DAA,.  The 
situation  is  somewhat  different  at  0  =  0,  as  seen  in  Figures  20  and  21.  In  Figure  20,  which 
shows  external-surface  pressure  histories,  the  DAA  humps  between  t  =  1  and  t  =  2  are  too  high, 
and  the  rise  times  immediately  after  t  =  2  are  too  long  relative  to  their  exact-history  counterparts. 
Even  so,  the  DAA2  history  captures  the  physics  rather  well,  exhibiting  accurate  peak  values  and 
correct  late-time  oscillation.  In  Figure  21,  which  shows  internal-surface  pressure  histories,  the 
DAA  histories  cannot  manage  the  deep  negative  excursion  of  the  exact  history  during  0.75  <  t 
<  2  and  the  abrupt  rise  at  t  =  2.  DAA2  performs  better  than  DAA„  producing  a  deeper  negative 
excursion  and  producing  correct  late-time  oscillation. 

Exact  and  DAA  shell-velocity  histories  are  shown  in  Figure  22  through  24.  DAA2  clearly 
outperforms  DAA„  the  latter  producing  histories  in  Figure  22  and  24  that  deviate  substantially 
from  the  exact  histories  at  late  time;  this  derives  from  DAA,’s  tendency  to  introduce  too  much 
acoustic  damping  (Geers,  1978).  DAA2  fails  to  reproduce  the  brief  oscillation  during  2  <  t  <  4 
in  Figure  22  and  23,  and  to  rise  rapidly  enough  at  t  =  2  in  Figure  24;  generally,  however,  it  does 
quite  well,  producing  good  general  response  with  accurate  peaks  and  correct  late-time  oscillation. 
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SECTION  9 

FIRST  ORDER  DAA  FOR  ELASTIC  MEDIA 


This  section  contains  formulations  of  first-order  doubly  asymptotic  approximations 
for  infinite  and  semi-infinite,  homogenous,  elastic  media.  The  formulation  for  the  infinite 
domain  is  a  generalization  of  the  matrix  form  given  in  Underwood  and  Geers,  1981,  and  then 
implemented  and  evaluated  by  Mathews  and  Geers,  1987.  Here,  the  method  of  operator 
matching  is  used  as  in  Sections  2  through  5  for  an  acoustic  fluid. 

First,  we  write  down  the  dynamic  extension  of  Somigliana’s  identity.  Next,  we  present 
the  first-order  early-time  approximation  for  infinite  and  semi-infinite  domains.  Then  we 
develop  the  first-order  late-time  approximation  for  the  infinite  domain,  which  is  followed 
by  a  similar  development  for  the  semi-infinite  domain.  Finally,  we  formulate  first-order 
DAAs  for  infinite  and  semi-infinite  domains  by  operator  matching. 

The  section  concludes  with  a  transformation  of  the  DAAs  from  operator  form  into  matrix 
form,  and  the  application  of  the  matrix  forms  to  some  canonical  problems.  Numerical  results 
are  compared  with  corresponding  “exact”  results  in  the  literature. 

9.1  DYNAMIC  SOMIGLIANA  IDENTITY. 

With  the  displacement  vector  u{x,y,z,t)  expressed  through  a  Helmholtz  decomposition 
in  terms  of  a  scalar  potential  <p{x,y,z,t)  and  a  vector  potential  tp(x,y,z,t)  as  (see,  e.g., 
Achenbach,  1973) 

u  -  V<p-\-V  x  ip,  V  ■  rp  =  0,  (9.1) 

the  wave  equations  for  a  uniform  elastic  medium  are 

c2DV7<p  =  if,, 


(9.2) 

c|vV  =  rp 


where  cp  and  cs  are  the  dilatational-  and  shear-  wave  speeds,  respectively,  given  by 


2  A  +  2fi 
CD  =  - - » 


(9.3) 
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An  exact  integral-equation  solution  to  (9.1)  and  (9.3)  for  the  displacement  field  on  the 
smooth  surface  or  surfaces  of  the  elastic  medium  is  provided  by  a  extension  of  Somigliana’s 
identity,  which  may  be  written  in  Laplace-transform  space  as  (Cruze  &  Rizzo,  1968) 


up (s)  +  J  uQ(a)fpQ{s)dSQ  =  J  tQ(s)UpQ(a)dSQ 


(9.4) 


In  this  equation,  P  and  Q  denote  points  on  the  surface  S,  up  and  tp  are  surface  displacement 
and  traction  vectors,  respectively,  and  Tpq(s)  and  Upq (,)  are  second-order  tensor  operators 
with  components 


Uii(s)  =  ^(#0  “  xR.iR  j) 


~  4ir[(rf£  RXXSii  dnR'jTli) 


(9.5) 


+<|  -^Tr-Tr-\*r^ 


where  the  usual  Cartesian  tensor  notation  applies,  where  the  PQ  subscript  for  Rpq,  defined 
after  (2.3),  has  been  omitted  for  notational  simplicity,  and  where  ip(R,s)  and  \{R,9)  are 
defined  for  three  dimensions  as: 


il>{R,s) 


(1  +  S  + 


cl  e-‘R/Cs 

s2R2  *  R 


Cjy 

32R2^  R 


(9.6) 


x{R,s) 


(i  + 


3cs 

sR 


+ 


34  _  c|  3cd  3 c*D 

s2R7  R  c2d{  sR  32  R2  R 


9.2  FIRST-ORDER  EARLY-TIME  APPROXIMATION:  ETA,. 

The  first-order  ETA  for  an  elastic  medium  is  a  vector  extension  of  that  for  an  acoustic 
medium  [see  (2.5)].  It  is  given  in  transform  space  for  an  infinite  domain  by  (Underwood 
and  Geers,  1981) 

ET A\  \  tp(a)  -  pC3Up(s)  (9.7) 
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where  C  is  a  second-order  tensor  that  involves  the  dilatational-  and  shear-wave  speeds 
defined  in  (9.3).  For  a  semi-infinite  domain,  (9.7)  also  applies  if  the  DAA  boundary  is 
somewhat  removed  from  the  flat  traction-free  surface  of  the  semi-infinite  domain,  or  if  it 
intersects  the  flat  surface  at  angles  equal  to  or  exceeding  90°;  we  restrict  ourselves  to  such 
situations. 

ETAi  is  a  local  approximation  ,  stating  that  each  element  of  the  surface  S  independently 
generates  three  plane  waves,  one  dilatational  and  two  shear,  which  propagate  normally 
outward  into  the  elastic  medium.  Because  it  approaches  exactness  only  as  s— ►  oo,  it  is 
singly  asymptotic. 


9  3  FIRST-ORDER  LATE-TIME  APPROXIMATION 
FOR  A  WHOLE-SPACE:  LTA’V. 

Here,  we  apply  to  (9.4)  the  procedure  described  in  Section  2.3,  keeping  terms  through 
s°  throughout.  This  yields  the  standard  Somigliana  identity. 


tip  +  J  UQTpQ{Q)dSQ  =  J  tQUpQ{0)dSQ, 


(9.8) 


where  the  tensor  components  are  given  by  [c.f.  (  9.4)  and  (  9.6)] 

1 


Uij(  0)  = 


47rp(l  -  v)R 


[(3  -  4 u)6ijR<iRj] 


(9.9) 


Z«(0)  =  8x(l  -  7)R?  {^((1  ~  2v)6ii  +  3R  ,Rj]  ~ (1  -  -  R  ^} 

The  static  relation  (9.8)  is  based  on  the  Green’s  function  for  an  infinite  elastic  medium 
obtained  by  Kelvin  in  1848.  With  the  spatial-operator  definitions 


BqQ  —  J  qQUpQdSQ 


(9.10) 


GqQ  =  JsqQ{6{P  -  Q)  +  TpqWSq 

where  6(P  -  Q)  is  the  Dirac  delta-function,  (9.8)  may  be  expressed  in  the  form 
LTA\V  :  tQ  =  B~xGtiP  (9.11) 

LTA]V  is  a  non-local,  singly  asymptotic,  quasi-static  approximation. 
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9.4  FIRST-ORDER  LATE-TIME  APPROXIMATION 
FOR  A  HALF-SPACE:  LTAf . 

The  first-order  LTA  for  a  semi-infinite  space  is  also  a  quasi-static  approximation  ex¬ 
pressed  by  (9.8),  but  with  Tpq  and  Upq  given  by  expressions  corresponding  to  the  half¬ 
space  Green’s  function  of  Mindlin,  1936,  rather  than  that  of  Kelvin.  These  expressions 
are  constructed  by  augmenting  (9.9)  with  additional  terms  (Brebbia,1984).  For  U ,  the 
augmentation  is 


UpQ  —  UpQ  +  UpQ 
in  which  the  components  of  Upq  are; 


U  n  =  Kd 
U\2  =  Kdr2 

U;x  =  Kdr2 
U22  =  Kd 


8(1  -  i/)2  -  (3  -  4v)  + 


(3  -  4i/)r2  -  2 cx  6cxrf 


R'2 


+ 


R'4 


(3  -  4i/)r1  4(1  -  i/)(l  -  2u)  6 cirj 


R'2 


R'  +  r[ 


R'4 


(3  -  4t/)ri  4(1  -  i/)(l  -  2v)  6 c£rj 

W2  +  R'  +  t[  ~  R'4  . 


r'2  2 cx  ( 1  -  3r^  \ 

~  +  ~R*  \  R'2  ) 

4r(l  -  u)(  1  -  2v)  /  _  r'2  \ 

*'  +  r'  \  R'{R'  +  r[)) 


II 

. *» 

Kdr2r3  [ 

3  = 

r3^12 

r 2 

Uii  = 

r*U21 

’"2 

U32  = 

^23 

^3*3  - 

Kd  1  + 

(3  -  4v)  4(1  -  i/)(l  -  2v)  6 cx 

R'2  R'{R'  +  r|)2 


(3  -  4v) 
R'2 


r'J  2 cx  (  3r?\ 

_  +  U2  \l  ~  ~Ri2j 


+ 


4r(l  -  i/)(l  -  2v) 
R'  +  r\ 


(-wfcl) 


(9.12) 


(9.13) 
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where  (see  Figure  25) 
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In  operator  notation,  LTAi  for  the  semi-infinite  domain  appears  the  same  as  that  for  the 
infinite  domain,  given  by  (9.11). 


9.5  FIRST-ORDER  DOUBLY  ASYMPTOTIC  APPROXIMATIONS  FOR 
WHOLE-  AND  HALF-SPACES:  DAAl 

Because  the  only  differences  between  the  first-order  early-  and  late-time  approxima¬ 
tions  for  infinite  and  semi-infinite  domains  reside  in  the  operators  Tpq  and  Upq,  we  can 
formally  develop  first-order  DAA’s  for  the  two  domains  simultaneously.  We  will  use  the 
method  of  operator  matching  for  this  purpose.  The  appropriate  trial  equation  is 

[sPi  +  P0]up(s)  =  tQ(s)  (9.18) 

where  the  spatial  operators  P0  and  P\  are  not  functions  of  s. 

For  s— *  0  we  write  (9.18)  as 

[P0  +  0(s)]uP(s)  =  tQ(s )  (9.19) 

and  match  with  (9.11)  as  s-+  0,  which  yields  Po  =  B~lG.  This  the  asymptotic  match  for 
the  static  limit.  For  s-+  oo  we  divide  (9.18)  through  by  s  to  get 

[Pi  +  0(s-1)]up(s)  =  s_1fQ(s)  (9.20) 

Now  we  match  with  (9.7)  as  s— *■  oo  to  give  Pi  =  pC.  Introducing  these  results  into  (9.18), 
we  obtain,  in  transform  space, 

DAAi(s) :  1q(s)  =  \pCs  +  R-1(5]  up(s)  (9-21) 

and  in  the  time  domain 

DAAi(t):  tQ{t)  =  pCuP(t)  + B~xGuP{t)  (9.22) 

Note  that  the  DAAi  for  elastic  domains  is  not  spatially  local,  because  of  the  late-time 
approximation  term. 
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9.6  MATRIX  DAAi  FOR  BOUNDARY  ELEMENT  ANALYSIS. 


The  most  direct  way  to  obtain  the  matrix  DAAi  for  either  whole-  or  half-spaces  is  to 
discretize  the  singly  asymptotic  approximations  (9.7)  and  (9.11  )  and  then  employ  the 
method  of  matrix  matching.  Thus,  we  preoperate  (9.11  )  through  by  B  and  introduce  into 
the  resulting  equation  and  into  (9.7)  the  finite  element  approximations 

tQ(t)  =  v£t(0 


UQ(t)  = 


(9.23) 


Then,  with  a  column  vector  of  weight-functions  wp  we  form  the  weighted-residual  equation? 


LTAi(f) : 


Bt(t)  =  Gu(<) 


(9.24) 


ETAi(t)  : 

t(t)  =  pCu(t) 

in  which 

B  =  j  vfpB\QdSp 

G  =  J  w-pGvqdSp  (9.25) 

C  =  J  WpCVqdSp 

If  we  now  follow  in  a  matrix  context  the  operator  -  based  matching  procedure  carried 

out  in  the  previous  section,  we  obtain  the  matrix  DAAi  in  an  elastic  whole-  or  half-space 

DAAi(t)  : 

t(t)  =  pCu(<)  +  B_1Gu(l).  (9.26) 

9.7  CANONICAL  PROBLEMS. 

It  is  useful  to  compare  DAA  based  and  “exact”  results  for  canonical  problems,  as  done 
previously  by  Underwood  and  Geers,  1981,  and  by  Mathews  and  Geers,  1987.  Here  we 
consider  two  canonical  problems,  both  pertaining  to  a  spherical  cavity  subjected  to  sudden 
internal  pressurization. 
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The  first  problem  is  a  spherical  cavity  embedded  in  an  infinite  elastic  medium  and  excited 
by  an  internal  step  pressure.  This  problem  possesses  radial  symmetry,  and  has  a  well-known 
analytical  solution  (Timoshenko  &  Goodier  1970).  With  a  as  the  cavity  radius  and  Pq  as 
the  pressure  magnitude,  the  radial  displacement  of  the  cavity  is  given  by 

u(t)  =  — jae~“‘(cos  ast  +  —  sin as't) 

4p  (,  s 

-ae~at(—  sin  as't  +  cos  as't)  (9.27) 

+  1  -  e_at(cos  as't  +  —t  sin  as't) 

where 

cD(  1  -  2u) 
a  a{l-v) 

(9.28) 

The  corresponding  analytical  DAA  solution  is  simply 

udaaW  =  ^  (l  -  .  (9.29) 

In  order  to  generate  numerical  DAA  solutions,  a  dynamic  boundary  element  program  has 
been  built  that  is  based  upon  the  program  constructed  by  Mathews  (Mathews  and  Geers, 
1987).  The  program  uses  eight  node  quadratic  quadrilaterals  for  spatial  discretization  and 
the  trapezoidal  rule  for  time  integration.  For  the  present  problem,  the  boundary-element 
model  for  the  cavity  boundary  consists  of  24  eight  node  elements  over  the  entire  spherical 
surface.  The  analytical  exact,  analytical  DAAi,  and  numerical  DAAi  solutions  are  shown 
in  Figure  26  for  the  parameters  p  =  1.00,  p  -  1/6,  v  =  1/4,  a  =  1,  and  Po  =  1.  The 
analytical  DAAi  and  numerical  DAAi  solutions  Eire  seen  to  be  almost  identical,  and  the 
DAAi  solutions  agree  well  with  the  analytical  exact  solution  at  both  early  and  late  times. 
As  previously  observed  by  Underwood  and  Geers,  1987,  the  DAAi  solutions  do  not  exhibit 
the  response  overshoot  seen  in  the  exact  solution.  This  is  characteristic  of  solutions  to  first- 
order  differential  equations  like  (9.26).  A  second-order  DAA  is  capable  of  accommodating 
such  overshoot. 

The  second  problem  is  a  spherical  cavity  embedded  in  a  semi-infinite  elastic  medium  and 
excited  by  an  internal  step  pressure.  Tliis  problem  does  not  posses  radial  symmetry,  and 
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does  not  posses  an  analytical  solution.  However,  a  boundary-element  solution  based  on 
numerical  inversion  of  Laplace  transforms  lias  been  generated  (Manolis  and  Ahmad,  1988), 
and  an  analytical  solution  to  the  related  static  problem  exists  (Bonefed,  1990).  In  terms  of 
the  geometry  shown  in  Figure  27,  the  latter  solution  is 


Ur  = 


-*>£-*£* 


i 

+  2 


uz  —  P0a' 


{2(1  - 


fl  -  — 

'  #-3  D.a 
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R13  R-3 


]} 


(9.30) 


where 


R'  -  yjr2  +  (2  -  z0)2 
R"  =  \Jr2  +  (z  +  zq)2. 


Numerical  DA  Aj  and  numerical  inversion  solutions  for  this  problem  are  shown  in  Figures 
28-30,  along  with  the  late-time  static  asymptotes  given  by  (9.30);  the  physical  parameters 
specified  are  the  same  as  those  previously  used  for  the  infinite-domain  problem.  Figure  28 
pertains  to  the  top  of  the  cavity,  i.e.,  the  point  on  the  cavity  surface  closest  to  the  free 
surface  of  the  elastic  half-space,  Figure  29  pertains  to  a  point  90°  around,  and  Figure  30 
pertains  to  the  bottom  of  the  cavity.  These  figures  show  that  the  DAAj  solutions  agree 
with  the  numerical  inversion  solutions  at  early  time  and  appear  to  approach  the  correct 
late-time  asymptotes;  unfortunately,  the  numerical  inversion  solutions  do  not  extend  far 
enough  in  time  to  allow  a  completely  satisfactory  comparison. 
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SECTION  10 


CONCLUSION 

This  report  documents  the  formulation  and  evaluation  of  new  doubly  asymptotic 

approximations  for  simplifying  the  analysis  of  transient  medium-structure  interaction  problems. 

More  specifically, 

1.  The  formulation  of  first-  and  second-order  DAA’s  for  an  external  acoustic  medium  has  been 
systematized;  finite-element  discretization  has  been  introduced  to  configure  the  operator-based 
formulation  for  boundary-element  solution. 

2.  First-  and  second-order  DAA’s  for  an  internal  acoustic  medium  have  been  systematically 
formulated  on  an  operator  basis;  finite-element  discretization  has  been  introduced  to  configure 
the  formulation  for  boundary-element  solution. 

3.  The  canonical  problem  of  a  spherical  shell  filled  with  an  acoustic  fluid,  submerged  in  an 
acoustic  medium,  and  excited  by  a  plane  step-wave  has  been  solved  by  modal  analysis;  special 
techniques  have  been  developed  and  applied  to  achieve  satisfactory  convergence. 

4.  Extensive  numerical  results  for  the  canonical  problem  have  been  generated  for  exact,  DAA, 
and  DAA2  treatments  of  the  internal  and  external  fluid-structure  interactions;  the  numerical 
results  have  been  compared  to  assess  DAA  accuracy. 

5.  The  formulation  of  the  first-order  DAA  for  an  infinite  elastic  medium  has  been  systematized; 
finite-element  discretization  has  been  implemented  to  configure  the  operator- based  formulation 
for  boundary-element  solution. 

6.  The  first-order  DAA  for  a  semi-infinite  elastic  medium  has  been  systematically  formulated 
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on  an  operator  basis;  finite-element  discretization  has  been  implemented  to  configure  the 
formulation  for  boundary-element  solution. 

7.  Boundary-element  DAA  results  for  suddenly  pressurized  spherical  cavities  embedded  in 
infinite  and  semi-infinite  elastic  media  have  been  generated;  the  DAA  results  have  been 
compared  with  corresponding  results  by  other  investigators. 

The  principal  conclusions  reached  during  this  study  are: 

1.  First-order  DAA’s  are  marginally  satisfactory  for  approximating  transient  medium-structure 
interactions  involving  external  and  internal  acoustic  domains  and  external  elastic  domains. 

2.  Second-order  DAA’s  are  highly  satisfactory  for  treating  external  acoustic  domains;  they  are 
satisfactory  for  treating  internal  acoustic  domains. 

3.  The  second-order  DAA  for  an  internal  acoustic  medium  is  sufficiently  accurate  to  warrant 
early  implementation  in  production  codes  for  underwater  shock  analysis. 

4.  Second-order  DAA’s  are  needed  for  treating  infinite  and  semi-infinite  elastic  media;  the 
techniques  used  herein  to  formulate  second-order  DAA’s  for  an  acoustic  medium  may  be  applied 
to  an  elastic  medium  as  well. 

5.  Further  DAA  development  is  desirable  in  order  to  obtain  approximations  of  higher  accuracy 
and  broader  application. 
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Figure  2.  Weighting  characteristics  of  modified  Ces&ro  summation  and 

standard  partial  summation  (CS3-N  =  Ces&ro  summation  over  modes 
3  through  N;  PSN  =  partial  summation  over  modes  0  through  N). 
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normalized  pressure  p/P, 


nondimensional  time  t 

Figure  4.  Incident-wave  pressure  histories  produced  by  modified 
Ces&ro  Summation  (CS3-N  =  Cesaro  summation  over 
modes  3  through  N,  e  =  m-s  error  over  0  <  t  £  2). 
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normalized  pressure  p/P, 


Figure  7.  External-  and  internal-surface  pressure  histories 
by  modified  Ces&ro  summation  (CS)  with  partial 
closure  (PC)  for  a  steel  shell  at  0  =  it 


Figure  8.  External-  and  internal-surface  pressure  histories 
by  CS  with  PC  for  a  steel  shell  at  0  =  0. 
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Figure  9.  Radial  shell-velocity  histories  by  CS  with  PC 
for  a  steel  shell  at  0  =  it  and  0  =  0. 


Figure  10.  External -surface  pressure  histories  at  0  =  it 
for  a  fluid-filled  shell  and  an  empty  shell. 
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Figure  13.  Radial  shell-velocity  histories  at  0  =  n 

for  a  fluid-filled  shell  and  an  empty  shell. 


Figure  14.  Radial  shell-velocity  histories  at  0  =  0 

for  a  fluid-filled  shell  and  an  empty  shell. 


9. 


Figure  15.  Pressure  and  fluid-particle-velocity 
histories  atr  =  0  inside  a  steel  shell. 


Figure  16.  Exact,  DAA,  and  DAA2  external-surface 
pressure  histories  at  6  =  n  for  a  steel  shell. 
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normalized  pressure  p/P(  normalized  pressure  p/Pi 


Figure  19.  Exact,  DAA,  and  DAA2  internal- 

surface  pressure  histories  at  0  =  vJ2. 


Figure  20.  Exact,  DAA,  and  DAA2  external- 
surface  pressure  histories  at  0  =  0. 
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Figure  25.  Three-dimensional  geometry  (in  the  case  of  the  half-space,  the  infinite  free 
surface  lies  in  the  X2  -  X3  plane). 
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Figure  26.  Radial  displacement  response  of  a  pressurized  cavity 
in  an  infinite  elastic  medium. 
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Figure  27.  Geometry  for  a  cavity  embedded  in  a  semi-infinite  elastic  medium. 


Figure  28.  Radial  displacement  response  of  a  pressurized  cavity 
in  a  semi-infinite  elastic  medium  (0  =  0°,  d  =  2a). 
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Figure  29.  Radial  displacement  response  of  a  pressurized  cavity 
in  a  semi-infinite  elastic  medium  (0  =  90°,  d  =  2a). 


Figure  30.  Radial  displacement  response  of  a  pressurized  cavity 
in  a  semi-infinite  elastic  medium  (0  =  180°,  d  =  2a). 
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UNIVERSITY  OF  COLORADO 
2  CYS  ATTN:  B  A  LEWIS 
2  CYS  ATTN:  P  ZHANG 
2  CYS  ATTN:  TL  GEERS 

WEIDLINGER  ASSOC,  INC 
ATTN:  H  LEVINE 

WEIDLINGER  ASSOCIATES.  INC 
ATTN:  TDEEVY 

WEIDLINGER  ASSOCIATES,  INC 
ATTN:  M  BARON 

WESTINGHOUSE  ELECTRIC  CORP 
ATTN:  D  BOLTON 
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